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Abstract. A new numerical method is introduced to study the problem of time 
evolution of generic non-linear dynamical systems in four-dimensional spacetimes. 
It is assumed that the time level surfaces are foliated by a one-parameter family 
of codimension two compact surfaces with no boundary and which are conformal 
to a Riemannian manifold ^ . The method is based on the use of a multipole 
expansion determined uniquely by the induced metric structure on c € . The approach 
is fully spectral — i.e. it avoids pointwise evaluations of the basic variables — in 
the angular directions. Instead, Gaunt coefficients as matrix elements are used to 
evaluate multilinear expressions. The dynamics in the complementary 1+1 Lorentzian 
spacetime is followed by making use of a fourth order finite differencing scheme. In 
handling the pertinent 1+1 transverse degrees of freedom the techniques of adaptive 
mesh refinement (AMR) is also applied. 

In checking the reliability and effectiveness of the introduced new method the 
evolution of a massless scalar field on a fixed Kerr spacetime is investigated. In 
particular, the angular distribution of the evolving field in superradiant scattering 
is studied. The primary aim was to check the validity of some of the recent 
arguments claiming that the Penrose process, or its field theoretical correspondence — 
superradiance — does play crucial role in jet formation in black hole spacetimes while 
matter accretes onto the central object. Our findings appear to be on contrary to these 
claims as the angular dependence of superradiant scattering of massless scalar fields 
does not show any preference of the axis of rotation. In addition, the characteristic 
properties of superradiance in case of a massless scalar field was also investigated. On 
contrary to the general expectations we found that by an incident wave packet, which 
had been tuned to be maximally superradiant, we cannot get more energy by the 
scattering process than the energy sent in. It was found that instead of the occurrence 
of energy extraction from black hole the to be superradiant part of the incident wave 
packet fail to reach the ergoregion rather it suffers a total reflection which appears to 
be a new phenomenon. 
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1. Introduction 

The basic equations of various theories are non-linear. In studying these types 
of dynamical systems analytic methods by themselves do not provide a completely 
satisfactory framework. Therefore it seems to be of fundamental importance to develop 
numerical methods that are capable to simulate long time evolution of non-linear 
dynamical systems. Motivated by this sort of necessities in general relativity various 
groups developed their codes aiming to make progress in the study of astrophysical 
systems containing black holes and neutron stars. Fully general relativistic simulations 
of coalescing binaries consisting of neutron stars and/or black holes are now possible by 
making use of variants of the generalized harmonic formulation [lj and moving puncture 
approach [2, 3j in the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism [U [3] 
(for a comprehensive review on the recent developments with additional references see, 

e.g., mm)- 

Besides these main stream efforts there are some apparently less ambitious ongoing 
projects trying to provide precise long term time evolution of various matter fields on 
fixed stationary background spacetimes which may or may not contain a black hole. In 
general, these investigations — due to the relative simplicity of the underlying physical 
system — provide an arena to test some of the new numerical methods before applying 
them to investigate the aforementioned much more complicated astrophysical systems. 
Immediate examples for these type of investigations with a stationary black hole as 
a background spacetime can be found, e.g. in [8J where high order finite differencing 
was applied, or in [9j where investigations of axial symmetric systems and the use of 
pseudospectral method (although only moderate angular momentum quantum numbers 
were involved), or in [TO] where results on the use of pseudospectral method without 
assuming axial symmetry but with utilizing parallel computing were reported. Similar 
dynamical systems were investigated in a series of papers [TTj [T2"l [T5| [HJ [TJ)|. In 
these papers the viability of the simultaneous use of the techniques of conformal 
compactification, along with the use of hyperboloidal initial value problem, in numerical 
simulations were demonstrated. 

In the present paper we introduce a numerical method to study the problem of time 
evolution of generic non-linear dynamical system in four-dimensional spacetimes. The 
time level surfaces are assumed to be foliated by a one-parameter family of codimension 
two surfaces which are conformal to a compact Riemannian manifold c it> without 
boundary. The degrees of freedom in directions tangential to 'rf — they are referred 
as angular directions — are treated with spectral representation (multipole expansion 
whenever c lo is homeomorph to a two-sphere S 2 ) which is based on L 2 expansions of 
the basic variables in terms of the eigenfunctions of the Laplace operator on c io . The 
fields in the transverse 1+1 dimensional spacetime directions are evolved by making use 
of the method of lines based on a fourth order finite difference numerical scheme. The 
pertinent numerical method incorporates the techniques of the adaptive mesh refinement 
(AMR). All the operations on the basic variables, involving angular degrees of freedom, 
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are done without applying point-wise evaluations, i.e. the method is fully spectral, not 
merely pseudospectral. 

There are various advantages ensured by this method. Firstly, the usual problems 
related to the coordinate singularities in the involved angular differential operators can 
be avoided. Secondly, all the operations in the angular directions which are linear in the 
basic field variables are exact. What is even more significant — by applying the Sobolev 
embedding theorem based arguments — all the non-linear operations such as pointwise 
multiplication or division by fields can also be treated within the spectral representation. 
The mathematical background of the applied new method — it is assumed to be known 
but seldom, if ever, collected in a systematic self contained way — is also presented in 
details in the appendices. 

In practice, all the multipole expansion series are truncated at certain finite order. 
In this respect the applied method is perturbative. Nevertheless, the error introduced by 
these approximations can be kept to be at a tolerable low level by increasing the number 
of the involved modes. The residual error is monitored and the precision, efficiency and 
the viability of the proposed method have been justified to be satisfactory. 

The introduced new method is applied to investigate the time evolution of a massless 
Klein-Gordon field on a fixed Kerr black hole spacetime. Within this setting, the angular 
dependence of the outgoing radiation was studied such that the initial data was fine 
tuned to have the highest possible potential for superradiance. Our investigations were 
motivated by some recent attempts trying to provide a physical model yielding high 
energy collimated matter streams (referred frequently as jets) originating from compact 
astrophysical objects. 

For instance, in |16j it was asserted that Penrose process involving Compton 
scattering on electrons and electron-positron pair production in photon-photon 
scattering can give rise to ejection of highly collimated matter streams along the 
axis of rotation. In [T7] an alternative support of these expectations was proposed. 
In particular, the possible existence of a class of timelike geodesies representing the 
worldline of the escaping particle yielded in the Penrose process — thereby emerging 
from the ergoregion — and having the axis of rotation of the black hole as an asymptote 
was examined. On contrary to these expectations in [18J where the evolution of a dust 
sphere falling onto a rotating black hole was considered no significant collimation effect 
had been found. 

In the present paper an analogous field theoretical model will be investigated. 
Distinguished attention will be paid to the angular dependence of outgoing radiation 
yielded by a scattering process with using initial data that has been fine tuned to possess 
the highest possible potential to generate superradiance. 

Let us also mention here that some preliminary studies of the dynamics of massless 
scalar field has been done by the present authors in |19| (see also |15j). However, 
neither the initial condition was fine tuned to generate to be superradiant solution, 
nor a detailed description of the applied numerical methods was given in either of 
these works. We would also like to emphasize that in parallel to the preparation this 
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paper the corresponding package of GridRipper with the implementation of the system 
investigated in this paper is made to be available for public use |20j. 

Let us also recall here that, based on the estimates in [2T1 [22] . superradiance 
is expected to be more significant whenever higher spin fields such as gravitational 
radiation is involved. Nevertheless, as stated above, throughout this paper 
considerations are restricted to the case of complex scalar fields. Similarly, our results 
concerning jet formation assume that the involved matter is modeled by a complex 
scalar field. Thereby our results do not exclude jet formation found in some recent 
astrophysically motivated more complex magnetohydrodynamical simulations (see e.g. 

[231 125). 

The paper is organized as follows. In Section [2j the physical setup including the 
field equations and coordinate choices are introduced. Section [3] to present an outline 
of the applied numerical method. Section [4] is to discuss some of the delicate issues 
related to the applied boundary conditions, while in Section [6] the initial data used 
in our numerical simulations is introduced. The main results are exposed in Section 
[7| while our concluding remarks are summarized in Section [8j The Appendices are 
to provide a systematic summary of the mathematical background of the applied new 
method, in particular, presenting all the details making it possible to use the techniques 
of multipole expansion in treatment of non-linear dynamical systems. 



2. Field Equations 

As mentioned above in this paper the evolution of a neutral massless scalar field 
propagating on the domain of outer communication of a fixed stationary Kerr black 
hole spacetime is considered. Although the code developed (which can be downloaded 
from [20]) is also capable to evolve a charged and self-interacting scalar field on a Kerr- 
Newman background in this paper attention will be restricted to the above mentioned 
simple case. We would like to mention that even this simple dynamical system is complex 
enough to test the viability and reliability of the proposed new method based on the 
spectral method. 

To start off let us recall first the Kerr metric given in Boyer-Lindquist coordinates 

t,r,$,tp [26J. The part of the spacetime on which our investigations will be carried out 

is the domain of outer communication that possesses the product structure IR 2 x § 2 and 

can be covered by Boyer-Lindquist coordinates t, r, $ and <p taking values from the 

intervals — oo < t < oo, < r < oo, < $ < it and < ip < 2tt. The metric g in these 

coordinates reads as 

A - a 2 sin 2 (tf) , 
g = — dt ® dt 

Zj 

a(r 2 + a 2 - A) sin 2 (i?) . , 
— (dt <8> dcp + dip <8> dt) 

Zj 

E T sin 2 ^) 

+ — dr <g> dr + £ d$ <g> d$ H — — dip <g> dip , gl) 

ZA Zj 
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where the smooth functions A, £ and T are determined by the relations 

A = r 2 + a 2 - 2 Mr, 
£ = r 2 + a 2 cos 2 ($), 

r = (r 2 + a 2 ) 2 -a 2 Asin 2 (tf). (§2) 

The symbols M and a denote the mass and the specific angular momentum of the Kerr 
black hole spacetime. The field equation of a complex valued scalar field <£> can be 
written as 

V a V a $ = 0, (2. 

which after the conventional first order reduction, for the vector variable $ r 
reads as 



d t $-- 
dt&t 



1 



(A 2 d r ^ r + 2 A (r — M) $ r + A L §2 ($) 
-a 2 d 2 ^-2a(r 2 + a 2 -A)d^ t 



d t <$> r = d r $ t 

where the differential operator 

1 



sin 1 !? 



d$ [sin ■& d$] + 



-—d 2 
sin 2 "d ^ 



(2 



5) 



is nothing but the Laplace operator on the unit sphere §> 2 with its canonical Riemann 
metric. To get rid of the coordinate singularity of the radial differential operator for the 
first multipole component of $ at the origin in the Minkowski limit — that can also be 
applied in other cases whenever an origin is present in the computational domain — the 
following conventional trick had been applied. Instead of $ the variable — r ■ $ was 
evolved using the field equation transformed accordingly. 

In the Kerr case with M > it turned to be rewarding to use instead of the r and 
ip the new ones r* and ip defined as 



r* r 



r + 



1 /ln(r — r+) ln(r — r 



+ 



K- 



In 



6) 



7) 



where 



1 r± 



2ri 



■± + a z 

is the surface gravity on the outer and inner event horizon located at r± 



(2 



(2 



(2 



M ± 



a 2 , respectively. By making use of these coordinates close to the event horizon 
much better resolution could be achieved which is supported by the fact that the null 
geodesies of minimal impact can be given as t ± r* = const, = const, (p = const. 
In consequence of the use of these new coordinates during the evaluation the inverse 
relation r = r(r*) had to be determined numerically which was done by implementing 
a simple Newton-Raphson method. 
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3. Numerical Evolution 

This Section is to provide a short outline of the applied numerical methods. As 
mentioned already the method is based on multipole expansion on each of the topological 
two-spheres determined by the t = const and r = const level surfaces. Accordingly, 
instead of evolving the fields $, $ t and $ rt themselves their multipole components [3>]J?\ 
[<&t\T an d [^rJr — which are functions of t and r* exclusively yielded by L 2 expansions 
of $, $j and $ rt with respect to the spherical harmonics — 

{Y?\e = 0,...,oo,m = -e,. ..,£}, §1) 

had been evolved. A comprehensive presentation of the mathematical background can 
be find in the Appendices, while the applied code GridRipper with the implementation 
of the investigated system can be found at |20| . 

3.1. Evolution in the t,r section 

The evolution equations for the multipole components [$]™, [<&t\T an d [^V*]™ were solved 
in the t — r* plane by making use of the 1+1 dimensional C++ based PDE solver of 
GridRipper described in details in [28, 29, 20j. The numerical algorithm utilized by this 
code is based on the method of lines in a fourth order Runge-Kutta scheme such that 
the spatial derivatives were evaluated with a fourth order symmetric finite difference 
stencil. To guarantee stability — by suppressing high frequency instabilities — a standard 
fifth order dissipation term, as proposed by Gustafsson et al [22] was also applied in 
solving the evolution equations for the multipole components [<&t]T an d [^V,]™- 

Note, that the use of this dissipation term does not affect the order of accuracy of the 
applied numerical scheme. 

The 1+1 algorithm of GridRipper makes use, as a built in package, the techniques 
of adaptive mesh refinement (AMR) as proposed by Berger-Oliger algorithm [30J (see 
also |29"t [28]). The use of AMR is based on the idea that a refining of the spacetime 
mesh has to be done at those locations where the Richardson error 

WfAtArjt^) ~ f2At,2Ar{t,r)\ 

2At (2« - 1) 

exceeds a predefined threshold, where fAt.Ar denotes the numerical solution obtained 
on a spacetime mesh with At temporal and Ar spatial finite difference, q is the order 
of accuracy of the finite difference scheme, and || • | is a semi-norm. As spatial and 
temporal refinement is performed simultaneously, the value of the Courant factor — i.e., 
the ratio of the temporal and spatial step size — remains intact, thereby in principle the 
stability of the finite difference scheme is not affected. 

In our simulations q took the value 4, while the semi-norm || • | was chosen to be the 
L 2 norm of the multipole expansion of c? r „\f f on each two-sphere. In order to be able to 
implement a relative error type quantity in specifying the tolerable error this L 2 norm 
of d r ^ was normalized by the pertinent L 2 norm of d r ^ on the initial slice. 
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3.2. Evolution in the angular if section 

The remaining angular (p directions were handled by a purely spectral method, 
completely avoiding point evaluation. The expansion coefficients [<&t]™ and [^r*]™ 

were stored in a C++ structure implementing an algebra defined by coefficient-wise 
linear operations and with pointwise multiplication of the basic variables. It is worth to 
be emphasized that the viable utility of the latter operation is not obvious at all as the 
multipole series, by construction, are guaranteed to be convergent only in the L 2 sense 
without an immediate support of their convergence in the pointwise sense. Therefore, 
in the generic case, the multipole expansion coefficients of pointwise products are not 
expected to be derived from the multipole coefficients of the factors without evaluating 
them pointwise and applying a subsequent numerical multipole expansion of the yielded 
product. Clearly, such a complicated approach would be computationally intensive not 
allowing the use of multipole expansions with sufficiently large £ values, e.g. £ > 16, to 
make the error introduced by truncation to be tolerably small. Nevertheless, whenever 
the basic variables are known to belong to the class of C 2 functions there exists a purely 
spectral approach that makes the evaluation of their pointwise multiplication possible. 
All of the underlying ideas — which are of fundamental importance in guaranteeing the 
effectiveness of the proposed new method — are justified with mathematical rigor in 
Appendix A and |Appendix B It is worth to be noted that every solution to a field 



equation involving second derivatives in the strong sense is of differentiability class 
C 2 . As the proposed use of the spectral method avoids pointwise evaluations there is 
a significant reduction in the required computational power in carrying out full 3+1 
dimensional simulations. This reduction is also supported by the fact that the Gaunt 
coefficients, introduced in Appendix B, which are necessary in evaluating products of 



multipole coefficients have to be calculated only once and stored them in the computer 
memory during the rest of the simulation. 

In practice, whenever pointwise products of C 2 fields truncated at l\ and £2 
multipole order is evaluated the result shall have non-vanishing coefficients up to 
l\ + £2 multipole order. Therefore, as opposed to linear operations, multiplications 
do not respect any prefixed maximal expansion order. Correspondingly, in the applied 
numerical approximation the multipole order of products have to be kept to be bounded 
which was done by truncating at the value max(£i, £2)- Note, however, that convergence 
tests have to be performed by varying the maximal allowed order £ max to justify the 
viability and the accuracy of the proposed numerical scheme. 

In evaluating the time derivative of the basic variables another critical non-linear 
operation has also to be performed. It is the division by a nowhere vanishing variable. 
The associated difficulties can be overcome by tracing back the operation of pointwise 
division to the operation of pointwise multiplication with the help of Neumann series 



expansions. This perturbative method as discussed in details in |Appendix C| and 



Appendix D| as it is another key ingredient of the proposed new method. It is also 



shown there that the necessary number of iterations in performing this perturbative 
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division method grows only with the logarithm of the required accuracy. 
3.3. Storage and computational requirements 

The storage of the basic variables via their multipole coefficients becomes even more 
efficient when the variables may be assumed to be of C°° class in the angular directions 
as in that case the corresponding sequence of multipole expansion coefficients are 
guaranteed to decay faster than any polynomial order as it is justified in |Appendix E 
The number of non-vanishing multipole coefficients of a variable, truncated at maximal 
order £ max , is (Cm + l) 2 - Therefore, the storage requirement is quadratic in £ max . For 
the number of the non-zero Gaunt coefficients the approximate formula 0.7 • (^ m ax) 4 ' 7 
can be verified numerically for £ max > 8. In particular, if the considered problem 
is axially symmetric the number of non-zero multipole coefficients is only £ max + 1, 
while the number of Gaunt coefficients necessary to evaluate non-linear terms scale as 
0.66 • (4nax) 2 ' 8 provided that £ max > 8. 

It is also informative to compare the computational expense estimates to that 
of alternative methods, for instance pseudospectral methods. These require pointwise 
evaluations of field values in angular meshes and subsequent evaluations of non-linear 
terms. In order to estimate the cost of such pointwise evaluation, one has to take into 
account that at each mesh point a sum over the indices £ and m has to be performed. For 
each mesh point in the $ angle coordinate a sum over the index £, with evaluations for 
the involved m values, consists of (£ max + l) 2 terms. In addition, the cost of a sum over 
of the m values — which may be evaluated by making use of fast Fourier transform — can 
be seen to go not better than 5(2£ max + 1) log 2 (2£ max + 1). Taking then into account 
that there exists £ max + 1 pieces of mesh points in as a minimal estimate of the 
total cost of a (non-approximate) pointwise evaluation in a pseudospectral method we 
get 5(£ max + l) 3 (2£ max + 1) log 2 (2£ max + 1). On the other hand, the evaluation of non- 
linear terms using matrix products in our method scales as ~ 0.7- (^, ax ). Therefore, the 
computational costs of these two methods appear to be comparable although in the range 
8 < 4nax < 32 there is about factor of ten preference on the side of the fully spectral 
method. Clearly, at the end of the simulation we also need to do pointwise evaluations 
in extracting the physical content of the yielded data. However, these evaluations need 
not to be done on each time level surfaces rather only at some specific ones. A slight 
additional advantage is that the involved vast number of matrix multiplications in the 
spectral case can be paralellized in a very effective way. Nevertheless, we admit that it 
is really the physical problem which should decide which method fits better. 

4. Boundary Conditions 

The proper treatment of the timelike part of the boundaries is of fundamental 
importance in both analytic and numerical evolutionary problems [3T| |32| |33| 134"] . On 
numerical side it is only a tiny technical part of the problem that the spatial derivatives 
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cannot be evaluated by making use a symmetric stencil close and at the boundary. The 
major part of the problems originates from the fact that there remains a freedom in 
specifying free data on the timelike part of the boundary [3T| 132] . 

In numerical approaches one of the most conventional treatment is to use the 
Sommerfeld outgoing radiation boundary condition. This method is based on the 
assumption that at the border the transformed variable v]> = r • $ has vanishing 
derivative along the outgoing radial null geodesies, which determines the value of \E^, 
at the boundary, in terms of the values of \1/ and ty r there. Note first that — given the 
simple form of the outgoing radial null geodesies in the Kerr spacetime in terms of the 
t, r* coordinates |27| — it is straightforward to implement this boundary condition in 
our numerical setup based on the spectral method. Nevertheless, in our first test runs 
the Sommerfeld boundary condition was found to yield unsatisfactory behavior in long 
term evolutions. One should keep in mind that even in the simple case of a scalar field 
on Minkowski background, in case of a non-spherical field configurations, only a much 
more sophisticated treatment [35J can provide a proper numerical treatment, and this 
approach does not generalize — at least not in a straightforward way — to more general 
background spacetimes such as the Kerr black hole. In virtue of this result one does 
not expect the Sommerfeld outgoing radiation boundary condition to work properly. 
Indeed, it was shown in that even in the simplest possible case of a massless 

Klein-Gordon field in Schwarzschild spacetime the asymptotic decay rate of the field 
may significantly be affected by the use of the Sommerfeld boundary condition. Our 
numerical experiments also justified (see Figure 4 1 below) that this outgoing radiation 
condition yields to significant instabilities at the boundary even in the short-term 
evolution of strongly non-spherically symmetric configurations. 

Another obvious idea is to carry out the numerical simulation near the boundary 
may be the following. Instead of applying any sort of outgoing radiation condition use 
the fourth order method of lines everywhere — as it is done in the interior — by making use 
of an asymmetric stencil close and at the boundaries. This simple minded approach also 
yields instabilities developing at the boundaries, although this occurs much later than 
in case of the Sommerfeld boundary condition. These type of instabilities most likely 
are consequences of the sum up of the error produced by the asymmetric fourth order 
stencils close and at the boundaries. While trying to cure this unfavorable behavior 
we invented the following simple trick. The order of the finite difference scheme was 
gradually decreased from 4 to 2 then to 1 such that we still had fourth order symmetric 
scheme at the last but two points, second order symmetric schemes at the last but 
one points and a first order upstream or downstream at the right or left boundaries, 
respectively. The first order asymmetric treatment at the boundaries could also be 
considered as a simultaneous combination of a linear extrapolation of the field variable 
next to the boundary point with the application of a second order symmetric differential 
scheme. Numerical experiments justified that this simple trick in evaluating the spatial 
derivatives close and at the boundaries — although with the price of a reduction of the 
numerical convergence rate there — stabilized the time evolution and, more importantly, 
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guaranteed a long term satisfaction of the energy balance relation even for the evolution 
of non-spherically symmetric field configurations. 



Figure 41 shows a comparison of the effectiveness of the Sommerfeld boundary 
condition, the simple 0(4) boundary condition — with fourth order asymmetric stencil 
close and at the boundary — and the developed 0(4 — 2 — 1) boundary condition — where 
the order of the finite difference scheme was gradually decreased from 4 to 2 then to 1— 
as described in details above. In particular, the time evolution of a rotating massless 
scalar field with a solid toroidal support is considered, by using either of these three 
boundary conditions, on Minkowski background. It can be seen that the energy outflow 
pattern at the outer boundary, located at r = 63, remains, in long term evolution, 
according to our expectations only for the case of 0(4 — 2 — 1) boundary condition. In 
the other two cases so much spurious energy flows back into the computational domain 
through the outer boundary that kills the evolution at t ~ 26 for the Sommerfeld 
boundary condition and at t ~ 100 for the 0(4) boundary condition. 



Boundary condition test 
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Figure [4}l. (Color online) The energy outflow as a function of temporal coordinate 
t at the outer boundary, located at r = 63, for the evolution of an initially rotating 
massless scalar field on Minkowski background is shown with the application of the 
Sommerfeld, 0(4) and 0(4 — 2 — 1) boundary conditions, respectively. It can be seen 
that in case of Sommerfeld or 0(4) boundary conditions a spurious energy flow back, 
from the outer boundary, occurs, while no such spurious energy flow back happens 
in case of the 0(4 — 2 — 1) boundary condition. These numerical investigations were 
carried out using £ max = 12 and 512 spatial points in the base grid such that n = 5 
AMR refinement levels were allowed. 

It is of obvious interest to know whether the 0(4 — 2 — 1) "boundary condition" 
proposed and used by us is a proper one. A sufficiently detailed investigation of this 
issue exceeds the frame of the present paper and the pertinent results will be published 
elsewhere. Nevertheless, a simpleminded explanation concerning the well-posedness of 
the associated initial-boundary-value problem will be given below. Before doing so let 
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us point to some of the most fundamental facts which should also support our claim that 
the proposed 0(4 — 2 — 1) differential scheme does impose proper boundary conditions. 

(1) While by applying the Sommerfeld condition the code crashes faster than in case of 
the 0(4) schema where no boundary condition at all had been applied. As opposed 
to this long term stability characterizes the use of 0(4 — 2 — 1) differential scheme. 

(2) The exponential convergence rate, along with the high precession of the energy 
and angular momentum balances (see Section [7]), could not be achieved without 
propagating all the physical modes towards the boundary such that they leave the 
computational domain without reflection, i.e., no spurious incoming modes are born 
at the boundaries. 

The simpleminded argument goes as follows: Assume that we have a first order 
system of hyperbolic field equation of the form u t = Au r + Bu for a vector valued 
field variable u. In applying the 0(4 — 2 — 1) scheme at the last grid point a first 
order downstream finite difference stencil is applied. This, however, can be seen to 
be equivalent to the application of a second order symmetric finite difference stencil 
combined with a linear extrapolation. These two operations guarantee that the second 
r-derivative, u rr , vanishes there. This, in virtue of this field equation, yields a mild 
restriction — of the type applied in Sections 9-11 of j25] — on the t-derivative, Ut, at the 
very last grid point, in spite of the fact that apparently only the field equations were 
imposed there. 

5. Treatment of an origin 

The presence of an origin on the time slices always requires a very careful and precise 
treatment. An origin shows up in various physically realistic situations. For instance, 
if the background is the Minkowski spacetime or in case of fully dynamical spacetimes 
containing a pulsating neutron star. It worth to be emphasized that the method outlined 
below is applicable not only in case of a Minkowski background but in the generic case 
of fully dynamical spacetimes, as well. 

Before proceeding note that in the fully dynamical situations the time slices might 
contain more than one origin or no origin at all (see, e.g. |2Z] for explicit examples). 
Nevertheless, in this Section attention will be restricted to the conventional single origin 
case. In order to avoid ambiguities the r-coordinate is arbitrary and the origin is assumed 
to be located at r = 0. 

In order to avoid the appearance of the usual r-coordinate singularity in evaluating 
the radial part of the Laplace operator at the origin, the new basic variable \I/ = r • $ 
has already been introduced. Whenever an origin is located at r = a symmetric 
fourth order stencil can be applied in determining the spatial derivatives at and in a 
neighborhood of the origin based on the observation that the multipole coefficients [\&]™ 
may formally be extended to negative radii according to the rule 




§1) 



Multipole Expansion in Time Evolution of Non-linear Dynamical Systems 



12 



This relation follows from the assumption that the original field variable $ is at 
least C 1 — thereby it is C 1 along arbitrary straight lines through the origin — and 
from the reflection property of the spherical harmonics Y™ under the transformation 
ip I — y 71 — ip -\- 7rJ§] 

In addition, if $ is assumed to be C 2 — this assumption should not be considered 



as extreme especially if one recalls that $ is subject to (2 3) — its spatial Laplacian 

.2) 



: <9 r 2 (r-$) + -^L s2 ($) (|5 



l -d 2 r {r •$) + 1 

has to be finite at the origin. From this, along with subsequent applications of l'Hospital 
rule, the relations 



m? = o, 

[$]f = o, d r [$}T = o f 



£ = 

e = i 

£>2 

follow. This, however, along with the substitution \1/ = r • $, implies that the relations 

£ = 0: [tf]™ = 0, = 0, 

£=1: [tf]? = 0, d r [^}f = 0, 

£>2: = 0, d r [*]? = 0, d 2 MT = (@4) 

hold for the multipole components ft?]™ at the origin. 

In consequence of the algebraic relations [^e] m = and <9 r [\I/]™ = 0, these hold for 
I 7^ at the origin, the field equations read there as 

dt[*]T = o, 

d t m? = o, 

dt[ * rli ~\ 0, otherwise, 

which are completely regular at r = 0. In spite of this apparently straightforward 
regularization of the singular terms at the origin a simpleminded numerical 



implementation of (5 5) still yield unstable evolutions for non-spherically symmetric 
configurations. A close look at the evolution justifies that numerical error starts to 
grow very quickly at the grid-point next to the origin. This can be understood by 
recalling that in the evaluation of the '0/0' type term 

1 

r 



„ 2 -M*) d 5 



6) 



§ Indeed, the radial derivative <9 r <I> of $(r, $, ip) corresponds to the directional derivative of $ along the 
radial direction characterized by the certain constant values of the angles d, ip. As $ is required to be 
C , the directional derivative d r & may, then, be numerically evaluated by making use of a symmetric 
stencil of the applied finite difference scheme such that values of <£> from both sides of the origin are used 
at and closed to the origin. Nevertheless, on the the opposite side the values of $ can be determined as 
*(r, TT-tf, ^+tt) = J2Zo YL=-l W Y?(v-#, ip+rr) = £~ Y^-t W (-l) £ Y?(ti, <p) . Thus, 
by formally extending the functions <J>™ to negative radii — by making use of the rule $™ H> (— 1)^$™ — 
the radial derivative <9 r $™ of the functions $™ can numerically be evaluated. 
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the higher multipole components acquire larger weight, which significantly magnifies 
the related numerical error. Nevertheless, this difficulty can also be overcome by a 
systematic application of the algebraic relations formulated by (51) and (54). By 
requiring these conditions to hold — where the first and second derivatives are assumed to 
be evaluated as dictated by the applied fourth order symmetric finite difference scheme — 
it turns out that the values of f^]™ at the origin and next to the origin are algebraically 
determined by the values of [^J™ next to next to the origin, with the only exception 



[\& r ]o which evolves according to the regular field equation <9t[\I/ r ]g = c^f^o [see (5 5) 
at r = 0. 



6. Initial Data and the Applied Grid 



A generic initial data specification to our evolution equations (2 4) is composed by three 
functions 0, <p t and r specified on the t = initial data hypersurface, denoted by So, 
such that beside the trivial constraint <p r = d r (j) for the corresponding solution $ the 
relations $|s = <f) an d $t |s = <fit also hold. It is straightforward to recast such an 
initial data specification for the rescaled field variable — r ■ $ which is given as a 
function of the coordinates t, r*, (p defined in Section [2j 



6.1. Superradiance 

Before proceeding and introducing our choice for the only freely specifiable functions 
ip and ip t on E let us recall some simple facts related to superradiance. First of all, 
as it was shown first by Carter in [38] in the coordinates t,r*,$,(p the dAlembert 
operator separates for the t-Fourier transformed field. More precisely, the temporal 



Fourier transform, of a solution $ to (2 4) may be decomposed as 

1 



vr 2 + a'' 



E E r w s ^m 



£=0 m= 



(6 



1) 



where co is the frequency in the time translation direction and S™ au) denotes the 
oblate spheroidal harmonic function with oblateness parameter au and with angular 
momentum quantum numbers £, m — they are eigenfunctions of a self-adjoint operator — , 
while for the radial functions R™ w a one-dimensional Schrodinger equation of the form 



drl 



+ 



ma 



to 



r 2 + a" 



A • VI 



Re 



(6 



with suitable real potentials V t 



can be derived from the field equation (2 4). 



Physical solutions to (6.2) are supposed to possess the asymptotic 



The conventional argument ending up with the phenomenon called superradiance 
goes as follows 
behavior 

-iuir, _|_ ^ g+iwr, 



R 



L.LU 



Te 



as r — )• oo 
as r — > r . 



(6 



3) 
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where Qh denotes the angular velocity of the black hole with respect to the 
asymptotically stationary observers [26], and with reflection and transmission 



coefficients, 1Z and T [SJ, respectively. Notice that this asymptotic behavior (6.3) 
presumes the existence of a transmitted wave submerging into the ergoregion. By 
evaluating the Wronskian of the corresponding fundamental solutions, "close" to infinity 
and "close" to the horizon, it can be shown that the reflection and transmission 
coefficients satisfy the relation [u — itlQh)T = (1 — TZ) oj [43] . Thereby, whenever 
1Z > 1 — or equivalently whenever the inequality < u> < mQn holds — positive energy 
is supposed to be acquired by the backscattered scalar wave due to its interaction with 
the Kerr black hole in the ergoregion. This phenomenon is referred as superradiant 
scattering which is also known as the field theoretical correspondence of the Penrose 
process derived in context of point particle mechanics |39| . 

6.2. Initial data 

In applying the introduced new numerical method our primary interest was to study the 
angular dependence of superradiant scattering. The applied initial data was specified 
accordingly — by applying an approach analogous to that of |4T)} 144] — and it was fine 
tuned to maximize the effect of superradiance. 

However, to investigate a clear manifestation of superradiance in a fully dynamical 
process — i.e. the way an incident scalar wave acquires extra energy by submerging into 
the ergoregion and then carrying it away from the central region — the initial data we 
applied is of compact support such that it is separated from the ergoregion on the 
initial time slice. Thereby, the initial data we have applied differs significantly, in its 
fundamental character, from that of [4THI41]. To fulfill the above mentioned requirements 
the initial data for the rescaled field variable ^ = r • $ was chosen as 

V(r,, Cp) = e-^~ r ^f{n - r* ) Y?{&, <p) , 

Mr*,&, 0) = <p) + e-'^-^fiu - no) *7(#, <f>) , 

Vv,(r*,tf,<£) = r ,V>(r„#,£) , @4) 

where / : R — > C is a smooth function of compact support, /' denotes its first derivative 
and Uo, r* are some real parameters. Note that the appearance of the extra r factor in $ 
has no effect on the above recalled argument concerning the appearance of superradiance. 
Indeed, this factor may be suppressed by redefining the function / that has not been 
specified yet. 

It can be seen that in an asymptotic region of the Kerr background (or everywhere 



if the background is the Minkowski spacetime) the initial data specification (6 4) yields 
an inward traveling spherical wave packet starting with a radial profile /(r* — r*o). 
Accordingly, in an asymptotic region the solution in a sufficiently small neighborhood 
of the initial data surface might be approximated as 

*(t, r„ tf, 0) w e-^ (r *- r * 0+t) /(r* - r* + t) y?(t?, <p). §5) 

It is informative to have a look at the temporal Fourier transform, J^J/, of this 
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approximate solution \P that reads as 

where a; being the temporal frequency while &f stands for the Fourier-transform of 
/. Assuming that & f — playing the role of a frequency profile function — is sufficiently 



narrow the approximate solution (6 5) looks almost like a monochromatic spherical wave 



solution similar in nature to the ingoing part of the wave determined by relations (6 1) 



and (6.3). Accordingly, by tuning Uq such that the energy flux absorbed by the black 
hole to become negative — this is expected to be achieved by choosing Uq such that 
< Uq < mVLn — one would expect that a to be superradiant solution is yielded. It can 
also be seen that the energy extraction may be maximized by choosing uq = hm£ln and, 
in addition, by guaranteeing that J™ H \JP f\ 2 (uj—uj ) dco « \ JPf\ 2 (uj— oo ) da;, which 
happens whenever the frequency spectrum is narrow enough to be entirely included by 
the superradiant frequency regime. 

We would like to emphasize that the above outlined construction of a to be 
superradiant initial data specification involves a number of heuristics assumptions. For 



instance, the Fourier spectrum (6 6) is assumed to represent the Fourier transform of 



the purely inward traveling wave (6 5) and whence the contribution from back scattering 



is completely neglected. To convince ourselves, in investigating the time evolutions of 
specific initial data choices, the power spectrum in temporal frequency of a supposed to 
be superradiant solution was also determined at a constant r* sphere which is located 
towards the black hole with respect to the compact support of the initial data. As it 



can be seen on Figure 7_6 the solution remains in the desired frequency regime. 

Based on the above outlined reasoning in our numerical simulations the radial profile 
function / :!-)•€ was chosen to possess the form 



ifx e [-- 



0. 



W Wl 

2 ' 2 J 

otherwise , 



©7) 



2 , 2 J. J-w 



which is a smooth function of the real variable x with compact support [— 
be compatible with our most important physical requirements that yields an incident 
wave packet that may acquire extra energy after penetrating through the ergoregion the 



initial parameter r * in (6 4) was chosen to be sufficiently large to have a clear separation 



of the support of the initial data and the ergoregion on S . 



6.3. Grid size and parameters 

The radial extent of the computational domain used in our simulations was chosen 
to be the closed interval —64 < r* < 64 in the massive case (M = 1), whereas the 
closed interval < r* < 64 in case of the Minkowski limit (M = 0). The specific 
angular momentum parameter a of the Kerr background was always chosen to be 0.9 
while the Schwarzschild limit was achieved by taking a = 0. The evolution of the 
system was investigated in the time interval < t < 192. The fine tuned parameters 
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of the initial data (6,4) — tuned to have the largest possible effect in superradiance — 
were co = 0.313394503136629, r* = 31.8229346475152, w = 35.3679317843828, while 
the angular and azimuthal mode numbers £ and m were fixed by choosing i = 2 and 
m = —2, 0, 2. Accordingly, the initial data had pure quadrupole character, while for m 
the values —2, and 2 signifies counter-rotating, non-rotating and co-rotating initial 
distributions, respectively. In virtue of the above discussion we may only expect the 
appearance of superradiance in the co-rotating case with m — 2, while no or negligible 
effect may be anticipated in the non-rotating or counter-rotating cases with m = or 
m = —2, respectively. 

In order to justify the above very specific choice made for the parameters oj Q , r*o 
and w let us recall the list of requirements they have to satisfy. 

• (u>/4) -1 <C rnQu ^ The width of frequency profile should be much smaller than 
the width of the superradiant frequency domain. 

• w <C r* )inax — r^ergosphere ^ The width of the initial data has to be much smaller 
than the part of the domain of outer communication outside to the ergosphere and 
covered by the grid. 

• w ^> Ar* The width and ramp of the wave packet has to be much larger than 
the spatial resolution of the base grid applied in AMR. 

• Uq <C At -1 f± The leading frequency of the initial data has to be much smaller 
than the maximal frequency allowed by the temporal resolution of the base grid. 

• r*o — \w > r* ergosphere ^ The support of the initial data has to be well separated 
from the ergoregion. 

• f*o + \ w < r *,max ^ The support of the initial data has to be included with suitable 
margins by the radial computational domain. 



6.4- Generic initial data for GridRipper 

Let us finally mention that in spite of the fact that in the investigations reported in this 
paper the initial data is always of pure multipole type in our code GridRipper (that can 
be downloaded from [20]) the generic case — whenever a multipole expansion of the initial 
data is required — is also implemented (see, e.g., |19j for an application). In the current 
version of GridRipper this is done by simply integrating numerically the product of the 
basic variables with Y™ over the r* = const angular spheres on S . In order to make 
this part computationally inexpensive — reducing thereby the required computational 
time to the order of seconds — the very efficient and precise two dimensional adaptive 
Genz-Malik (AGM) method (36) is applied. 



7. Numerical investigations 

This Section is to introduce our main results concerning the evolution of a massless 
scalar field on Kerr background. As emphasized earlier distinguished attention will 
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be paid to the angular dependence of the field and to the formation of superradiance. 
Before presenting our numerical results it is important to justify the reliability of the 
proposed new method. 



7.1. Error estimates and convergence 

As emphasized in Section [3] the representation of the basic variables by truncated 
multipole series can only be 'exact' in the case of linear field equations. On the 
other hand, whenever the evolution equations contain non-linear expressions of the 
basic variables — with non-trivial angular dependencies — the multipole method becomes 
inherently perturbative. Nevertheless, it is believed that the error yielded by the 
truncation of the infinite multipole series remains at a tolerable level provided that 
the value of £ max is kept at a sufficiently high value. In order to demonstrate that 
this expectation is valid, the £ max dependence of some estimates on the error and the 
convergence will be shown below. 

Almost all of our simulations were performed by using £ max = 12, nevertheless, in 
order to be able to determine the convergence rate simulations with £ max = 14, 16 and 
18 were also performed in the Kerr case with initially co-rotating and counter rotating 
distributions. Note that as the initial data had pure quadrupole character the indicated 
variation of the value £ max had no effect on it. 

In what follows the numerical representation — with maximal multipole order £ max — 
of a function / will be denoted by f iia ■ Assume that A£ max is some positive integer. 
As a measure of the relative error of the variable /, the quantity 



II f - f II 



max. ■ - i ' max 



1} 



was applied. Notice that E (f) monitors the time dependence of the difference 

of the basic and finer solutions f, and f. , relative to the finer one. The norm 



|| ■ || applied here, and in (7 2) below, is the C° norm bounded from above by the second 
Sobolev norm C||- H/j-a^c), where C is the minimal Sobolev constant associated with the 
if|(S 2 ,C) C C°(§ 2 ,C) Sobolev embedding as discussed in Appendix A and Appendix 



pi Clearly, E (f) <C 1 has to hold for reasonable numerical solutions provided 

that the value of £ max is sufficiently large. 

Another useful quantity characterizing the validity of the applied numerical schema 
is the d'Alembert convergence factor defined for the numerical representation f t as 



II f - f II 

Q I' J Imax + Almax ^ lmax + 2Al max II 



' Ui;o,~^' )n:i.\ 1 1 f /' 



max + ^^max ' 



2) 



max- 



This quantity measures the local convergence rate centered at £ max + A£ n 

In virtue of d'Alembert's criterion guaranteeing a sequence to be summable, 

convergence of the numerical solution in £ max occurs provided that the inequality 
limsup, ^ (ft mj!) A<™(/)) < 1 holds. 
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A related quantity — which is useful in quantifying the appropriateness of the 
numerical scheme — is the local exponent of convergence defined by the ratio 



In (Q t 



if) 



(7 



3) 



"'-max 

As discussed in Appendix E[ whenever a function / is C°° the convergence has to be 



faster than any polynomial in £ max , and therefore the local exponent of convergence 
becomes an informative measure of convergence rate — which is expected to be constant 
whenever the convergence is exponential. 

On Figure 7 1 the time and £ max dependencies of the relative error E 
relevant for the basic variable \l/ are shown for an initially co- or counter rotating massless 
scalar field on a Kerr background with parameters M = 1 and a = 0.99. As it is clearly 
visible the relative error is always smaller then 10 -8 and it is decreasing while £ max 
is increased. If considerations are restricted to the initial part of the evolution, i.e., 
to the part before the wave packets leave the computational domain at t ~ 75, the 
relative error does not exceed the level ~ 10~ 12 . Following this initial, truly dynamical, 
period eight order smaller amplitude quasi-normal ringing and finally an even smaller 
amplitude power low tail decay occur (see, e.g., [IS]). The amplitude of these processes 
is comparable to the applied accuracy of the present simulations which yields a visible 
increase in the relative error. Notice also that the apparent linear hierarchy of the 
graphs — transparent in the applied logarithmic scale in both of the subregions — justifies 
that, indeed, the convergence rate in £ max is exponential. 

Let us mention here that as the finite differencing part of our code, applied in the 
"t — r" Lorentzian sector, is exactly the one which, along with the AMR part, went 
through careful and detailed convergence tests — and the pertinent results can be found 
in |28l [29] — we would like to recall here only that this part of the code is of fourth order 
accurate as it should be according to the implemented numerical scheme. 

The computational times as listed in Table |7|1| — relevant for the same systems 



as described in connection with Figure 7 1 and also for the common PC architecture 



AMD Phenom(tm) 2.3GHz CPU — justify that the required computational resources are 
affordable, i.e. the proposed new method is computationally inexpensive. 



7.2. Energy and angular momentum balances 

In addition to the rate of convergence in £ max the use of the energy and angular 
momentum balance relations provides another important consistency check verifying 
the reliability of the proposed numerical algorithm. 

Recall that the balance laws relate values of energy and angular momentum on 
portions of t = const hypersurfaces to energy and angular momentum fluxes across the 
timelike hypersurfaces connecting the edges of them. In particular, the argument goes 
as follows. Whenever there is a divergence free vector field J a on a spacetime it can be 
justified by referring to Stokes' theorem that for a spacetime domain N with boundary 
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t 

Figure [TJl. (Color online) The coordinate time and £ max dependence of the relative 
error E of "3/ representing the evolution of initially co- or counter rotating 

massless scalar field with £ = 2, to = 2 or £ = 2, to = —2, respectively, on a 
Kerr background with parameters M = 1, a = 0.99 and for the particular values 
f max = 12, 14 and 16 are shown. The apparent linear shifting of the error curves on 
this figure verify that the rate of convergence in £ nULX is exponential. 



9 

<-max 


Comp. time (m = —2) 


Comp. time (m = 2) 


12 


39582 sec 


56374 sec 


14 


53274 sec 


76089 sec 


16 


68900 sec 


99248 sec 


18 


87514 sec 


124727 sec 



Table [7}l. The computation time of the evolution of initially co- or counter 
rotating massless scalar field with £ — 2, to = 2 or £ — 2, m = —2, respectively, 
on a Kerr background with parameters M — 1, a — 0.99 within the coordinate 
time interval < t < 192 and with the particular choices of the values 
f max = 12, 14, 16 and 18. The indicated times were measured by using PCs 
with architecture AMD Phenom(tm) 2.3GHz CPU which justify that the proposed new 
method is computationally inexpensive. 



dN and outward pointing unit normal vector n a at ON the balance relation 

n a J a = [ V a J a = (7.4) 

I dN JintN ]_ 

holds. On the other hand, it is well-known that — as the vector fields dt and are 
Killing vectors on a Kerr spacetime — the contractions = —T a j,d^ and J£ = T a bd^, 
which are the energy and angular momentum currents, are divergence free, where Tab 
denotes the energy-momentum tensor of the matter fields. In our investigations, N(t) 
was chosen to possess — in the tortoise Boyer-Lindquist coordinates — the form of the 
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Figure [7j2. (Color online) The time dependence of the relative variation SE = 
■gjj- Jq n m n aJp and SL — Jgjy(t) n a^£ °f energy and angular momentum balances 
during the evolution of a massless scalar field on Minkowski, Schwarzschild and Kerr 
background spacetimes with initially co-rotating (m = 2), non-rotating (m = 0) and 
counter-rotating (m = — 2) pure quadrupole type initial data. The reference values Eq 
and Lo are the initial energy and angular momentum contents of the selected parts of 
the initial data surface t = 0, respectively. 



The constant r* values determining the edges of spatial section of the cylindrical 
domain of integration N(t) — in order to keep some margin from the edges of the 
computational domain — were chosen to be such that r*i = —63 and r 3f2 = 63 for the 
Kerr or Schwarzschild cases with M — 1, whereas r*i = and r*2 = 63 were used in the 
Minkowski limit with M = 0, where r* reduces to r. 



On Figure |7|2| the time dependence of the relative variation SE and 5L of energy 
and angular momentum balance relations are shown. Here SE and SL are defined as 

If 1 

SE = — n a J% and SL = — 
Eo JdN Lo 

where Eq and Lq denote the energy and angular momentum of the initial configuration 



n aJL ) 
dN 



(7 



5) 



within the spatial region r* G [r*i,r*2]- The graphs on Figure 7.2 make it transparent 
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that the energy and angular momentum balances hold up to a remarkable precision for 
the entire evolution. 

Note, finally, that the evaluation of the involved integrals can be done in a 
straightforward way in context of the spectral method as the integration with respect 
to the angular degrees of freedom can simply be given as L 2 scalar products of the basic 
variables which can be evaluated as outlined in Appendices A, B, C and D. In addition, 
the integrals with respect to the radial and temporal directions were evaluated by using 
a fifth order integration scheme to avoid the loss of accuracy of the numerical data 
yielded by the applied fourth order finite difference scheme in the t — r* plane. 



7. 3. Angular dependencies of the fields 

After presenting the consistency checks of the applied numerical scheme let us turn to 
the description of the physical properties of the solutions. In this Section our main 
concern is the angular dependence of the evolving scalar field. 

To have some hints regarding the dynamics of a massless scalar field with co-rotating 
quadrupole type data on Kerr spacetime the energy density and the momentum current 



distributions are shown on Figure |7|3| on the initial data surface, at t — 0, and on 
an intermediate time level surface after a scattering of the inward falling radiation has 
happened, at t = 48. It is visible that, in spite of the fact that the initial data was also 
fine tuned to be maximally superradiant, the dominant part of the outgoing radiation 
leaves the central region without indicating the slightest preference of directions close 
to the axis of rotation. 



On Figure 7 4 the time dependence of both the total integrated fluxes of the radiated 
energy and angular momentum through the r* = 63 sphere and the fluxes of the 
radiated energy and angular momentum integrated on the caps of sphere yielded by 
the intersection of a double right circular rotationally symmetric cone with apex angle 
29 = 7r/3 and the r* = 63 sphere is shown for various configurations. 

It can be seen that in the case of rotating initial configuration (with m ^ 0) the 
outgoing radiation is suppressed in the vicinity of the axis of rotation. An effective 
evaluation of the flux integrals on the caps of the sphere yielded by the intersection of 
a double right circular rotationally symmetric cone in the spectral framework requires 



additional technicalities which are described in details in |Appendix F 



It still remains really cumbersome to extract some insight concerning the anisotropy 



of the outgoing radiation simply by inspecting plots of the type depicted by Figs. 7 3 



and 7 4 Therefore it is important to have a clear measure of anisotropy. Assume that 
J a is a conserved current and consider a ball ^(f*) of radius r* = f* in the outer region 
of the computational domain. Denote by 33(f*,6) the disjoint union of the two caps 
yielded by the intersection of a double right circular rotationally symmetric cone with 
apex angle 29 and the ball 3${f*). Clearly then 3${f*) = SSif*, 7r/2). Finally, denote 
by A"^ ut | e (t) the integral S\oi\x3S{f e) nT * a ^ a ^ wnere n r,a stands for the unit form field 
normal to the hypersurface [0, t] x ^(f*, 8), i.e. n Tta points to the increasing r* direction. 
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Figure [t|3. (Color online) The spatial distribution of energy density n^ b T a b d^ and 
the r,d,Lp components of the energy current T a b d^ are shown at t = (top panel) 
and at t = 48 (bottom panel) for a massless initially co-rotating scalar field on a 
Kerr background with parameters M = 1 and a = 0.99. Note that only the sections 
corresponding to the azimuthal slices ip = and 7r are plotted. The initial data 
was fine tuned to be maximally superradiant and by t — 48 a scattering has already 
happened. The energy density, n^i 1 T a b d t a , is indicated by the color map while the 
drbT a b d^ and d^f,T a b 9 t a components of the energy current tangent to the ip = and it 
plane are indicated by arrows, whereas the azimuthal component dipbT b df is depicted 
by isocurves. The location of the singularity, the event horizon and the ergosphere is 
also indicated on central parts of the plots. Note that for the sake of simplicity the 
quantities indicated are given by referring to the Boyer-Lindquist coordinates. 
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Figure [7j4. (Color online) The time dependence of the integrated energy and 
angular momentum fluxes during the evolution of an initially quadrupole type co- 
rotating (m = 2), non-rotating (m = 0) and counter rotating (m = —2) massless scalar 
field on Minkowski, Schwarzschild and Kerr spacetimes is shown. Besides the total 
integrated fluxes the fluxes integrated on the disjoint caps of the sphere with radius 
r. t = 63 yielded by the intersection of a double right circular rotationally symmetric 
cone with apex angle 26 = ir/3 and the sphere are plotted. Note that in the Kerr case 
the co-rotating initial data with m — 2 was fine tuned to be maximally superradiant. 



Based on the above introduced quantities, as a measure of anisotropy, we may use 
then the expression 

4tt ^outlTr/2 l^maxj 

which is nothing else but the ration of the time dependent angle average of the integrated 
flux of the current J a through the two caps, ^(f*, 9), of the ball S3{f*) of radius r* = f* 
located at the north and south poles and of the angle average of the total integrated 
flux through the entire ball ^(f*). 
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Notice that, for any fixed 9 G (0,7i"/2) value, whenever the radiation has no 
anisotropy at all £/ J (t,9) tends to 1 as t — > i max , it tends to a value smaller than 1 
if the radiation shows preferences of the directions close to the equatorial plane, while 
£^ J (t, 9) tends to a value greater than 1 if the radiation prefers the axial directions. 
Clearly by choosing 9 to be small the sharp preference of the axis can be tested, whereas 
by increasing its value the anisotropy can be tested for a wider range of directions around 
the axis. In all of our investigation 9 was chosen to be 7r/6, i.e. all directions within a 
right circular rotationally symmetric cone with apex angle 29 = 7r/3 were included. 

On Figure [7 _5 the time dependence of the energy and angular momentum radiation 
anisotropics, £/ jE (t,7t/6) and L (t, 7r/6) is shown for Kerr, Schwarzschild and 
Minkowski background spacetimes. The evolution starts with co-rotating or counter 
rotating initially field configurations. As before the co-rotating initial data was fine 
tuned to be maximally superradiant. By the inspection of Figure 7_5 the following 
simple observation can be made. 

• Regardless whether the initial configuration was co-rotating or counter rotating a 
strong preference of the directions close to the equatorial plane is justified by the 
asymptotic behavior of E (t, 7r/6) and £/ jL (t, vr/6) both of which tend to a value 
much smaller then 1. 

• It is also clearly visible that the properties of the background spacetimes have no 
noticeable effect on the evolution of either £/ jE (t,7r/G) or £/ jL (t,7r/6). 

• Finally, the fact that the initial data for the co-rotating configuration was fine tuned 
to have the solution to be superradiant has no effect at all on the anisotropy of the 
energy and angular momentum distribution of the outgoing radiation. 



7.4- Total reflection and superradiance 

Based on the observed insensitivity of the anisotropy on the superradiant or non- 
superradiant character of the initial configuration it turned to be important to 
understand the reflection and absorption processes during the evolution of the 
investigated scalar field on black hole backgrounds. 



As it was already indicated in subsection |6.2| it is important to be sure that the 
type of initial data chosen there does correspond to be superradiant configurations. To 
justify that this is indeed the case we determined the temporal frequency spectrum of a 
numerical solution with initial data parameters u>q = |mfi#, r*o = 31.823, £ = 2, m = 2 
and with Kerr background parameters M — 1, a — 0.99 at the location r* = 14 which is 
located towards the black hole with respect to the compact support of the initial data. 
It is clearly justified by Figure [7[6] that the spectrum is indeed well contained within the 
superradiant regime as expected. 



The generic behavior of the incident wave packets is depicted on Figure 7 7 This 
figure shows the time dependence of the radial coordinate densities of the energy and 
angular distributions of the massless scalar field evolving on Minkowski, Schwarzschild 
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Figure[7j5. (Color online) The time dependence of the energy and angular momentum 
radiation anisotropics. s^ Je (t, ir/6) and g/' ,L (t, ir/6), is shown for Kerr, with M = 1, 
a = 0.99, Schwarzschild, with M — 1, a — 0, and Minkowski, with M = 0, 
a = 0, backgrounds. The evolution starts with co-rotating or counter rotating initial 
field configurations. The co-rotating initial data was fine tuned to be maximally 
superradiant in the Kerr case. 



and Kerr background spacetimes. The initial data is of quadrupole type and co-rotating 
or counter rotating in the Kerr case, according to the choices £ = 2 and m = ±2. 

The radial coordinate density of energy and angular momentum are the quantities 
$ and Jzf with the help of which the energy and angular momentum, E and L, on a 
t = const time level surface can be given as E = f t=const £dn and L = J t=conat ^dr^, 
i.e. in $ and Jzf the energy and angular momentum densities are integrated with respect 
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Figure [7j6. (Color online) The power spectrum in temporal frequency, at the location 
r* = 14, of an inward traveling wave packet with I = 2, m = 2 and with initial data 
having leading frequency u>o — hmilff. Thereby, the incident wave packet is, as it 
is expected, to be maximally superradiant. The parameters of the Kerr background 
were M = 1, a = 0.99. Note that the spectrum appears to be relatively intact in its 
superradiant character, i.e. the solution remains in the desired frequency regime. 



to the angular degrees of freedom, and they also involve the not yet integrated part of 
the 3-volume form induced on the t = const time level surfaces. 

We have found that in case of a massive background with M > for non- 
superradiant type of initial configurations, as it is expected, considerable part of the 
incident wave packet gets to be absorbed by the black hole. However, for initial data 
fine tuned to generate a totally superradiant configuration — contrary to the generic 
expectations — no energy extraction from the black hole was observed. Instead a total 
reflection, as it is shown on the bottom right panel of Figure [7/7 occurred. Notice 
the similarities characterizing the evolution of the to be superradiant initial data in the 
Kerr case (lower right panel) and the evolution of the scalar field on simple Minkowski 
background (top left panel) with no black hole in the setup. 

In trying to figure out the significance of total reflection recall now that to be able 
to produce extra energy in the process of superradiance part of the radiation — after 
submerging into the ergoregion — has to descent towards the event horizon. The lack of 
energy extraction from the black hole can now be understood as the reflection happened 
before the radiation could have reached the ergoregion. 

As our result concerning superradiance is on contrary to the conventional 
expectations it is important to emphasize that the observed phenomenon of total 
reflection, for the part of the wave packet belonging to the superradiant regime, was 
found to be robust with respect to the variation of the parameters of the background 
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Figure [7j7. (Color online) The time dependence of the radial coordinate density of 
energy and angular momentum, S and S£ (for their definition see the main text), is 
shown for Minkowski, Schwarzschild and Kerr background spacetimes. The evolution 
starts with co-rotating or counter rotating initially field configurations. As before the 
co-rotating initial data was fine tuned to be maximally superradiant in the Kerr case. 



spacetime and that of the initial data. Note also that as our pertinent results appear 
to be inconsistent with the claims [40, 41 J it is important to clear up the reason beyond 
these controversial conclusions. In doing so start by recalling that there is a significant 
difference between the type of initial data applied in |4TJ| |4"T] and in this paper. While 
the initial configuration we applied is of compact support in |4TJ| |4"T] the initial data 
was arranged to have non-trivial values everywhere in the ergoregion. Moreover, it is 
claimed in j4TJ| |4"T] that energy extraction from the black hole does occur. However, in 
our checks — applying horizon penetrating slices as in [TH] and exactly the same type of 
initial data as in |4TJ1 |4"T] — the to be superradiant character of the field was lost in an 
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extremely short period. More importantly, no energy flux leaving the ergoregion could 
be observed. All of these observations should be completed by emphasizing that our 
conclusions are not at all incompatible with claims in |4~2 l Pf3l l4~4"] . First of all, although 
in jl3] compactly supported data is applied in deriving analytic estimates concerning 
superradiance the yielded results therein are converted, on page 833, to quantitative 
estimates based on approximations derived by Starobinskii |32], in spite of the fact that 
the approximations applied in [42J are not entirely compatible with the use of compactly 
supported initial data. What is even more important is that the pertinent conclusion 
in jl3] provides only an upper bound for the gained energy which is about ~ 1% for 
the case I = m = 2. It may sound strange but this upper bound may be considered 
to be entirely compatible with the total reflection observed in our simulations as 0% is 
definitely smaller than this rough upper bound. 

In providing some more convincing evidences let us emphasize first that the 
conventional arguments of Misner and Zel'dovich supporting the existence of energy 
extraction are based on the use of individual modes. Note, however, that the study of 
the linear stability problem for Kerr spacetimes |44j . with the application of finite energy 
wave packets, taught us the lesson that statements at the level of individual modes need 
not imply statements for the superposition of infinitely many modes. 

Figure 7 8 is to provide an additional justification of our main result. On this 
figure the time dependence of the energy extraction coefficient, S °~J? out is shown, where 
S and S out , respectively, stand for the initial value of the energy and for the integrated 
energy flux through the ball of radius r* = 63 located at the outer boundary of the 
computational domain. As it is visible the graph of S °~f o ° ut starts at the value one and 
tend to zero from above. Note that the slowly decreasing part, with 150 < t < 1000, 
represents only the beginning of a long lasting quasi-normal ringing of the scalar field on 
the black hole background. The energy stored in these ringing modes will eternally be 
also radiated to infinity (see, e.g., [E]). Since S °~J? out does not change sign the energy 
radiated to infinity remain always smaller than So, which justify our conclusion that no 
energy extraction had happened. 

Clearly, one could claim that any numerical method has its own limitation which 
is true also in the present case. Nevertheless, in virtue of Figure 7 2 the accuracy of 
our numerical scheme allows us to put sharp upper bound on energy extraction which is 
~ 10~ 4 times So that is significantly smaller than the ~ 1% of So derived by analyzing 
individual modes |2Tj |22l |12] . 

Let us finally mention that the power spectrum in temporal frequency of the solution 
provides some new insight what happens whenever the to be superradiant wave packet 
approaches the ergosurface. On Figure 7 ^9] the r* dependence of the power spectrum in 
temporal frequency of the solution, which had been averaged for the angular degrees of 
freedom, is shown. The initial data is exactly the same quadrupole type with i = m = 2 



as used to generate the solution depicted on Figure 7.7 (Note that to determine 



the proper Fourier transform of the solution for all the indicated values of r* it was 
necessary to evolve the initial data both forward and backward in time.) It is clearly 
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Figure [7J8. (Color online) The time dependence of the energy extraction coefficient, 
^°^f out ^ s snown f° r quadrupole type initial data, where Sq and S ont , respectively, 
stand for the initial value and for the integrated energy flux through the sphere of 
radius = 63 located outward with respect to the support of the initial data. 



visible that the solution stays in the superradiant regime not only in the distant region 
but up to the ergoregion. In addition, it is also important that the frequency of the 
involved modes grows up to reaching the value u = m Qh where, in virtue of the relation 
(u — itlVLh) T = (1 — 1Z) oj, inevitably a total reflection has to occur. 

Let us close this Section by commenting the non- negligible reflection visible on the 
top right panel of Figure 7 7 depicting the evolution of the massless scalar field on a 
Schwarzschild background. One might be surprised by this reflection as intuitively it is 
tempting to assume that the Schwarzschild black hole would be ready to absorb almost 
the entire of the incident wave packet. Recall, however, that the angular momentum 
balance relation does not support the occurrence of such an overwhelming absorption. 
In addition it is also worth to have a look at the equation on a Schwarzschild background 
governing the evolution of a pure Y™ mode. In fact, the t — r* part of the wave equation 
for the coefficient ty™ = ty™(t,r*) reads as 

{%-%,)*? + Vi(r)*? = 

with the potential Vi(r) = (l — ^y-] (^—^- + ~rl> The repulsion — responsible for the 
reflection of the inward falling radiation and, in turn, leading to the celebrated power 
decay law of Price |Tf] H8] — is transparent on Figure 7 10 depicting the potential Vi(r) 
for the I = 0, 1,2 cases. What is really important here is that the maximum value of 
the potential V^5^ ~ 0.24 is significantly larger than ~ 0.096 which, in virtue of the 
types of arguments contained e.g. by Section II of Chapter III in [35], can be used to 
justify the observed scale of reflection. 



Multipole Expansion in Time Evolution of Non-linear Dynamical Systems 

Killing-time power spectrum of solution 



1 1 \ 


i i i 

M = 1, a = 0.99, m = 2 










Superradiant regime ' 
End of ergosphere: r = 2M ' 

i i : 


1 1 1 



10" 



10 5 



10 4 



'9. 
-a 

XI 



-64 



-39 



-14 



11 



36 



61 



10 3 



30 



Figure [7^9. (Color online) The dependence of the power spectrum in temporal 
frequency of the solution, which had been averaged for the angular degrees of freedom, 
is shown. The applied initial data is exactly the same quadrupole type with £ — m = 2 
as used to generate the solution depicted on Figure |7|7| The solution stay in the 
superradiant regime for the entire evolution while the frequency of the involved modes 
grows up to reaching the value uj = mQu where total reflection has to occur. 



8. Summary 

Our main concern in writing up this paper were at least two folded. On the one hand, 
we intended to introduce the generic setup of a method that is expected to provide a 
powerful new tool in studying the problem of time evolution of non-linear dynamical 
systems in four-dimensional spacetimes. On the other hand, we applied the introduced 
new method to study the evolution of a specific dynamical system. More precisely, the 
evolution of a massless scalar field on a fixed Kerr spacetime was investigated such that 
distinguished attention was paid to the angular distribution of the evolving field and to 
the occurrence of superradiance. 

In spite of the fact that the mathematical background of the introduce new 
method — which is mainly contained by the appendices — makes it to be applicable to 
dynamical systems the time level surfaces of which can be foliated by a one-parameter 
family of codimension two surfaces which are conformal to a compact Riemannian 
manifold c (o without boundary, in most of the cases with time level surfaces possessing 
the topology of R 3 or R x § 2 it suffices to assume that ^ is topological two-sphere § 2 . 

One of the main advantages of the new method is that it is fully spectral — 
not merely pseudospectral — in the angular directions. Thereby, whenever the basic 
variables are guaranteed to be at least of class C 2 in the angular directions, the spectral 
components — even though non-linear expressions of the basic variables are involved — 
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Figure [7}l0. (Color online) The potential Ve(r) involved in the wave equation over 
Schwarzschild spacetime is plotted for values i = 0,1,2. The potentials have their 
common zero value at the horizon, r = 2M, while they attain their maximum close to 
but on the domain of outer communication side of the event horizon. The occurrence 
of the partial reflection on a Schwarzschild black hole background can be understood 
by taking into account the repulsing character of this potential, along with the fact 
that the maximum value of the potential Vj^f ~ 0.24 is significantly larger than 
ujI « 0.096. 



can be evolved without their steady pointwise evaluations. Accordingly, the angular 
degrees of freedom — directions tangential to c io — are treated by applying I? expansions 
of the basic variables in terms of the eigenfunctions of the Laplace operator on c io . The 
corresponding expansion coefficients of a basic variable are evolved in the transverse 
1+1 dimensional spacetime directions by making use of the method of lines based on a 
fourth order finite difference numerical scheme such that the adaptive mesh refinement 
(AMR) is also incorporated. 

The main advantages associated with the use of the proposed new method are: 

• the coordinate singularities associated with the angular differential operators are 
treated in a fully analytic way 

• all the non-linear operations such as multiplication of the basic variables or the 
division of an expression by a nowhere vanishing variable — the latter can traced 
back to multiplication with the help of the Neumann series expansion — can be 
treated within the spectral representation without steady pointwise evaluation of 
these expressions 

• a very effective treatment of an origin of the time level surfaces can also be done 
by making use of fully analytic considerations. 

A systematic self contained presentation of the mathematical background of the applied 
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results and the implemented elements of the spectral method can be found in the 
appendices. 

The introduced new numerical method is used to study the time evolution of a 
massless Klein-Gordon field on a fixed Kerr black hole spacetime. We would like to 
emphasize that this dynamical system is already complex enough to explore the main 
technical elements of the proposed numerical framework. In particular, two of the most 
important ones, i.e. the evaluation of various multilinear expressions and the division 
by a function based on the use of spectral method were applied in carrying out all of 
our simulations. 

As all the multipole expansion series are truncated at certain finite order, £ max , the 
proposed numerical method is perturbative. By studying the error introduced by the 
involved approximations it was verified that the error can be kept at a tolerable low 
level by applying sufficiently many members of the multipole expansions. In addition, 
a suitable notion of convergence was also introduced, based on d'Alembert's criterion 
guaranteeing the summability of sequences. It was justified that the convergence is 
exponential in the value of £ max . 

In studying the time evolution of a massless Klein-Gordon field on a fixed Kerr 
black hole distinguished attention was paid to the precise characterization of the angular 
dependence of the outgoing radiation, as well as, to the development of superradiance. 
Our main related results are as follows. 

• There are attempts (see, e.g. [16, 17J) aiming to provide a simple and viable physical 
explanations of high energy collimated matter streams originating from compact 
astrophysical objects by applying Penrose process, or superradiance. On contrary 
to the underlying speculations we have found that the outgoing radiation has no 
preference at all of the axial directions regardless whether the initial configuration 
was to be superradiant or a more generic type. 

• In studying superradiance — in the particular case of a massless scalar field — we 
intended to investigate its formation by using an incident wave packet, which was 
fine tuned to be strongly superradiant, from the outer region of the domain of outer 
communication. On contrary to the general expectations, we found that instead of 
energy extraction from the black hole the incident wave packet failed to reach the 
ergoregion. Thereby, instead of superradiance, according to the authors knowledge, 
as a new phenomenon a total reflection occurs before the "to be" superradiant part 
of the incident wave packet reaches the ergoregion. A complete characterization of 
the new phenomenon definitely exceeds the scopes of the present paper. We plan to 
carry out its further investigations, involving possibly more generic type of initial 
data and, more importantly, a mixture of analytic and numerical techniques. 

Concerning the required computational resources it is also important to be 
emphasized that various implementations of pseudospectral methods even in case of 
the study of time evolution of a massless Klein-Gordon field on a fixed Kerr black hole 
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require powerful computers and parallel computing. As opposed to these the proposed 
new method was found to be very effective in the sense that time evolutions of the 
very same dynamical systems could be done within reasonable computational time on 
a stand-alone average personal computer without making use of parallel computing. 
We would also like to mention that the GPU implementation of the proposed method, 
which is under development, promises a further significant boost in the optimal use of 
the computational resources which will be inevitable in studying the evolution of more 
complicated non-linear dynamical systems. 

Finally, we would like to emphasize that our code GridRipper, together with the 
implementation of the of time evolution of a massless Klein-Gordon field on a fixed 
Kerr background and other examples, is made available for public use and it can be 
downloaded from |2Uj . 
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Appendix A. Theory of Multipole Expansion of Multilinear Expressions 

Let us start by recalling some of the basic notions and notations that we shall 
apply throughout the succeeding appendices. By a Riemannian manifold M we shall 
always mean an n-dimensional paracompact manifold (possibly with boundary) of 
differentiability class C r endowed with Riemannian metric g, which will usually be 
suppressed. We shall denote by L 2 (M, C) the Hilbert space of square-integrable complex 
valued functions^]] on a Riemannian manifold M, and by C l (M,C) the vector space of 
/-times continuously differentiable complex valued functions on M, where < I < r. 
Similarly, C l b (M, C) stands for the Banach space of bounded C l (M, C) functions on M 
equipped with the C l supremum norm, while (M, C) denotes the Banach subspace 
of C l b (M, C) functions on M which have zero limit at infinity!^ 




|| Here in a precise formulation instead of 'functions' function equivalence classes should be used, where 
two Lebesgue measurable functions are considered to be equivalent if they are almost everywhere equal. 
Referring to these function equivalence classes as simply 'functions' is a common practice in functional 
analysis which we also shall use. 

% A C l b (M,C) function is said to have zero limit at the infinity, if for any monotonously growing 
sequence of compact sets covering M, the C l supremum norm of / over the complement of the compact 
sets tends to zero. 
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Consider now the vector space of r-times weakly differentiable complex valued 
functions over M. The Sobolev norm of a function / belonging to this linear space 
is defined as 



£ / |V«/| 2 , (A.l) 
1=0 Jm 



were | • | denotes the pointwise norm generated by the Riemann metric g, V denotes the 
Levi-Civita covariant derivation determined by g, while the associated volume form is 
suppressed in the applied notation. It is important to keep in mind that the highest, 
r th -order, derivation in the above formula is required to be defined in the weak sense. 
The subset of the vector space of r-times weakly differentiable complex valued functions 
over M comprised by elements whose Sobolev norm || • ||_h-2( M)C ) is finite is called the 
Sobolev space and it is denoted by H%(M, C). The Sobolev space together with its norm 
forms a Hilbert space, since its norm is generated by an inner product, and also it is 
complete with respect to this norm. 

Let n, r and k be non-negative integers, such that r > | + k. The classical result 
of Sobolev embedding theorem (see e.g. |19]) asserts then that the relation 

H?(R n ,C) c C^(R n ,C) (A.2) 

holds, and also that there exists a positive real constant C such that the inequality 

II " ||C*,(R™,C) < C || • ||if2( R n )C ) (A. 3) 

is satisfied. These assertions are also known to hold |49j if W 1 is replaced by a compact 
subset in IR n . In the next part of this appendix, our aim is to provide a simple and 
self-contained justification of the fact that the Sobolev embedding theorem may also be 
applied in case of compact Riemann manifolds]^ For yet another alternative reasoning 
see, e.g., Ref. [50] . 

In doing so consider first a finite dimensional real vector bundle, W(M), of 
differentiability class C r over a manifold M that belongs to the same differentiability 
class. Then, a pointwise mapping of C r -sections of W(M) onto C°-sections of M x R 
will be called to be a C r -norm field if its pointwise restrictions to the fibers give rise to 
norms Indeed, every finite dimensional C r - vector bundle over a paracompact manifold 
M admits C r -norm fields. To see this recall that a norm field can always be defined 
locally over a coordinate chart, e.g., by taking in every point the natural Euclidean 
norm defined by a trivialization. These locally defined norm fields may, then, be sewn 
together by making use of a partition of unity subordinate to a locally finite collection 
of coordinate charts on M. The norm fields over a vector bundle are equivalent as it is 
justified by the following lemma. 

+ Note, however, that there are non-compact Riemann manifolds such that the Sobolev embedding 
theorem does not apply to them. 

* A more adequate way of formulating this definition is that a W(M) — > M x R C° fiber bundle 
homomorphism is a C r norm field if its restrictions to the fibers are norms. 



Multipole Expansion in Time Evolution of Non-linear Dynamical Systems 



35 



Lemma 1. If | • | and | ■ |' are C r norm fields over W(M), then there exists a positive 
real C r field C, such that \ ■ | < C\ • |. 

Proof. The proof is based on the paracompactness of M and on the equivalence of norms 
on a finite dimensional vector space. 

To start of consider a locally finite atlas = {(Ui, tpi) \ i G /} of M with partition 
of unity {Fi\i G /} such that each Ui has compact closure, denoted by Ui, in M. 
Assume that the dimension of the fibers of W(M) is N. Let us fix a trivialization 
(ey | j G {1, ... , N}) of W(M) over each particular chart (Ui, (fi) G s$ ' . 

As a consequence of the equivalence of norms on a finite dimensional vector space, 
for any p G M there exists a positive number c p , such that | • | < c p | • | p . Furthermore, 

c p may be chosen to be sup j^p. 

s p &W p (M)\{O p } |Splp 

By making use of the trivialization (e^j) of iy(M) over (U,ifi) it can be verified 
immediately that 



sup sup -j — p J = sup 



( N ' \ 

I E S 3 e i,j\ (p) 

sup 

v SeRjV ' |s|=1 iE^.|(p) 



(A.4) 



holds for any choice of (Ui,(fi) G The right hand side of (A.4) is a finite positive 
number, because it is nothing but the maximum of a positive valued continuous function 
over the compact manifold Ui x S^" 1 , where S N_1 denotes the A^ — 1 dimensional unit 
sphere. Let us denote this positive number by q. Then, 



<a\-\ (A.5) 



holds over [/,-. 



As an immediate consequence of ( A.5 ) we have that J~i\-\ < CjJ-i | • | holds throughout 
M, in accordance with the fact that J 7 ; is non-negative and supp(J-i) C Z7j. Note, then, 
that the sum Eie/ i s a positive valued C r function — which, as a consequence of the 
local finiteness of srf ' , has only finite non-zero terms in a sufficiently small neighborhood 
of any point in M — , and that by definition Eie/^i = 1- These, along with the above 
observations, implies then that | • | < (Eie/ c i^i) I ' I holds, which justifies the assertion 
of the lemma. □ 

We shall also apply the following two lemmas in verifying that the Sobolev 
embedding theorem may also be applied in case of compact Riemann manifolds. 

Lemma 2. Let (\ ■ |/)/e{o...., m } be norm fields, and V, V' be two C ' -covariant derivative 
operators. Then, there exists a positive C '-function C over M such that 



£iV' ( °-|*<C>r|V«.|<. (A.6) 

1=0 1=0 
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Proof. The assertion of this lemma may be justified by combining the following sequence 
of simple observations. 

(i) the covariant derivation V' can always be expressed as a sum of terms involving V 
and the C r_1 class Christoffel symbols, 

(ii) the triangle inequality holds for norms, 

(iii) the composition of a norm with a linear map is a semi-norm, 

(iv) the sum of a norm and a semi-norm is a norm, 

(v) and, finally, by taking into account Lemma [7J 

□ 

Now, as a direct consequence of Lemmas [T] and [2] we have the following. 

Lemma 3. Let (| ■ |z)ze{o,...,m} and (| • |;)jefo,...,m} be norm field collections, furthermore, 
V and V' be two covariant derivative operators as above. Then, there exists a positive 
C r function C on M such that 

j:\V' il) .\;<C±\V^.\, (A.7) 

1=0 1=0 

Being armed with the above results we can turn to the generalization of Sobolev 
embedding theorem to compact Riemann manifolds. 

Theorem 4. The Sobolev embedding theorem applies to compact Riemann manifolds, 
i.e., whenever M is a compact Riemann manifold of dimension n and r > | + k for 
non-negative integers r and k, the relation 

# r 2 (M,C)cCt(M,C) (A.8) 

holds, and also there exists a positive real constant C such that the inequality 

II • He* (m,c) < C II • \\h?(m,c) (A. 9) 

is satisfied. 



Proof. Obviously, the inclusion (A.8) holds as an immediate consequence of the relevant 
differentiability assumptions, and because C£,(M,C) = C 6 fc (M,C) = C k (M,C) due to 
the compactness of M. 

To justify ( A.9[ ) choose an arbitrary finite atlas srf = {(Ui,ipi)\i G 1} of M with 
partition of unity {Fi\i G /} subordinate to it. In proceeding, choose a real number e 
such that < e < 1, and denote by Vf the pre-image of the interval [e, 1] by the map J-i 
for each i e I. Then, by applying the conventional Sobolev theorem over each compact 
set Vf C supp(J-i) C Ui we have that 



sup 



^ \i=o J i=o Jv ^ 




J> W /I <a£/ |V«/| 2 (A.IO) 



for any / G H^(M,C), where now | • |, V and the volume form are assumed to be 
determined by the Euclidean metric associated with the local coordinates, <fi(Ui) C M n , 
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on Ui. In virtue of Lemma^ one may replace the Euclidean | ■ |, V and the volume form 
with the norm, covariant derivation and the volume form determined by the Riemannian 
metric g on M with the understanding that, as a compensation, the values of Sobolev 
constants Ci have to be adjusted accordingly. Furthermore, as Ti > e over Vf, we also 
have that 

2 



sup 

Vf 



£>w/|) W ^|v«/i 2 , (A.ii) 

v 1=0 J i=o Jv i 

where the monotonicity of the integral of a non-negative function has been taken into 
account. By making use of the above observations we also have that 



£ / iv w /i 2 = EE / *iv w /i 

J M ;r-T ,_n JUi 




2 



^Ea su ? Ei v(!) /i • ( A - 12 ) 

iei % v i \i=o J 
This is exactly the point where we utilize the compactness of M. Accordingly, for the 
rest of the proof we shall assume that the index set / is finite, which immediately implies 
that for some positive constant C e 

holds. By choosing e such that < e < 4r we also have that ^jVf = M since otherwise 
there would be a point p G M so that J-"i(p) < e for all i 6 /, This, however, is 
impossible since then ^ie/^O ) < Sie/^— ^ wou ld hold on contrary to the definition 



of the partition of unity. This, in virtue of (A. 12) and (A.13), justify then that (A.9) 



holds. □ 

In applying the above results consider now a compact Riemann manifold M. As 
it is well-known [50] the eigenfunctions of the Laplace operator are in C£(M, C), and, 
more importantly, their linear span is dense in any of the Banach spaces C l b (M, C), 
with < I < r. In addition, C l b (M, C) C Hf(M, C) is dense and — in consequence of the 
Holder's inequality — the C l supremum norm is stronger than the Hf norm over compact 
manifolds. Therefore, the linear span of the Laplace eigenfunctions is dense in Hf(M, C), 
and, in particular, in L 2 (M, C). It is also known that a linearly independent eigensystem 
of the Laplace operator over a compact manifold may be chosen to be orthonormal with 
respect to the L 2 scalar product, hence these form a complete orthonormal system in 
L 2 (M, C), a complete orthogonal system in H 2 (M,C), as well as, a Schauder basis in 
C l b (M,C), with < / < r. Let / G H 2 (MX) be some function and {Y l } ieX be an L 2 - 
orthonormal eigensystem of the Laplace operator. Then, as the elements of the system 



Multipole Expansion in Time Evolution of Non-linear Dynamical Systems 



38 



{Yi}i& are eigenvectors of the Laplace operator, in virtue of the Gauss theorem, we 
have that 

to f) H ?(M, C ) = (ibi-^y) to ( a - 14 ) 

for each i G X, where Aj denotes the associated eigenvalue, and DM = has also been 
assumed. In order to simplify some of the succeeding expressions we introduce the 
function S r : C — > C as 



(-Z) r+1 -1 ■€ , _-, 

S r (z) = Y(-z) 1 = I ' 2 ^ ' (A.15) 

1 r + 1, otherwise. 



1=0 



As an immediate consequence of (A. 14) the vector system {Yj,/ ^/S r ( Aj) is 



orthonormal in H 2 (M,C). Therefore, for the series expansion of the function / in 
H 2 (M,C) with respect to the complete orthonormal system {Y^/ a/ SV(Aj)}j S 2 

e(t='/) • 7 = = Eto/Wo^ (A.i6) 

is satisfied, which relation implies then that the L 2 series expansion of / with respect 
to {Yi} i£ x is also convergent in the H 2 sense. 

By combining the above observations we have that whenever the assumptions of 
the Sobolev embedding theorem holds then the multipole series expansions — which are 
convergent in the L 2 sense — will also be convergent in the uniform C k sense. 

Corollary 5. Assume that f is a function that belongs to C r (M,C), where r > | + k. 
Then, the multipole series expansion of f with respect to an eigensystem {Yi} ie x of 
the Laplace operator is also convergent in the C k sense. Furthermore, the convergence 
is independent of the summation order, as the system {Yi} i£ x shall remain to be a 
complete orthogonal system in L 2 (M, C) and in H 2 (M, C) after any index permutation. 
Therefore, pointwise absolute convergence of the multipole series also follows for the 
derivatives up to the order k, as the absolute convergence in a finite dimensional Banach 
space (which is nothing but C in the present case) is equivalent to summation order 
independent convergence. 

Remark 6. It is worth keeping in mind that the Schauder basis property of the Laplace 
eigenfunctions in C k (M, C) does not guarantee that an arbitrary function / e C k (M, C) 
can always be expanded in the form of / = ^2 ieX fiY%- (If that was true the relation 
fi = (Yi, f) L 2, i € X, would immediately follow from Lebesgue's theorem of dominated 
convergence.) Our statement follows from the fact that whenever {Y^i G X} is a 
Schauder basis that guarantees merely that the linear span of this system is dense 
in C k (M, C), i.e. every field / may be approximated by a sequence of finite linear 
combinations of this basis. This means that the series expandability with respect to 
a Schauder basis follows automatically only in Hilbert spaces, but not in general in 
Banach spaces. 



Multipole Expansion in Time Evolution of Non-linear Dynamical Systems 



39 



It is important to be emphasized that the pointwise absolute convergence property 
plays a crucial role in non-linear problems. To see this note that whenever the conditions 
of the above corollary are guaranteed to hold the pointwise product / • g of multipole 
expansions / = ^2 ieX jiXi and g = Yliexdi^i may be written in the optimal form 



f-9 




YY 



(A.17) 



= f* 9 J 

where the sums are understood in the pointwise manner, which, due to Fubini's theorem, 
are absolute convergent along with their derivatives up to the order k. The sums may 
also be understood in the uniform C k manner independently of the summation order, 
however, in this case we do not get absolute uniform C k convergence. 

The most significant advantages associated with the use of multipole analysis 
manifest themselves in evaluating non-linear terms. In fact, whenever the pointwise 
absolute convergence is guaranteed for / • g we have that 



f-9 = 52(f-g) h Yk 



(A.18) 



kex 



where 




2 fa 

iex jex 



YkYiYj . 



(A.19) 



M 



In deriving (A.19) we have used that, in virtue of Lebesgue's dominated convergence 



theorem, the summation and the integration can be interchanged, and also that the 
double sums are known to be uniformly convergent. Therefore, whenever the matrix 
elements jYiYjY^ with k G X, which are also known as Gaunt coefficients, can 
be determined, then, the multipole series expansion on multi-linear expressions may 
constructively be evaluated by making use of the multipole expansion of its factors. 



Appendix B. The Explicit Form of Gaunt Coefficients on Two-sphere and 
on n-dimensional Torus 

It is known (see, e.g., |51| [52] ) that the Gaunt coefficients on a two-sphere, S 2 , may 
be given either in terms of the Wigner 3j symbols or in terms of the Clebsch-Gordan 
coefficients. It is also known that the Clebsch-Gordan coefficients can be evaluated by 
various numerical algorithms, e.g., by the Racah formula (53] (see also |54|). 

In providing the Gaunt coefficients let us denote the spherical harmonics, which 
comprises the familiar orthonormal basis, by Y™ with integers £, m satisfying the 
relations < i < oo and — £ < m < I. 

Then, the Gaunt coefficient J g2 YJ^Y^Y 1 ^ 3 is known to be zero if either of the 
following conditions holds 
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(i) m 1 +m 2 + m 3 ^ 0, or 

(ii) 4< |4 -4j | or 4 +£ 2 < 4, or 

(iii) 4 + 4 + 4 is odd. 

The Gaunt coefficients may also be related to the Wigner 3j symbols as 

/ v mi V m2 V m 3 

Js 2 

/ (24 + i)(24 + i)(24 + i) (444W4 4 4\ fB ) 

V ~ 4tt " ' Vo o oy v^i^2m 3 y ' 

where the Wigner 3j symbols are given by the Clebsch-Gordan coefficients as 

3 = 1 • (4 4 mxma |444 (~m 3 )). (B.2) 

\mi m 2 m 3 / ^24 + 1 

The Gaunt coefficients relevant for the case of an n-dimensional torus, T n , are 
known to possess an even simpler structure. By making use of the orthonormal 
eigenstates Yfci,...^™, labeled by the integers k 1 , . . . , k n of the Laplace operator on T n , the 
integrals J Jn Y k \ ^n -Y k \ ^.n -Y k i^ k n may be seen to be zero if there exists % G {1, . . . , n} 
such that k\ + k\ + k\ ^ while it is J- otherwise. 

Appendix C. The Division as an Operation within the Framework of 
Multipole Expansions 

One of the most delicate issues while using multipole expansions arises whenever we 
have expressions containing the operation of division by a function. Clearly, if this 
function vanishes somewhere, then, a direct evaluation of the division is not tolerated 
by any of the numerical methods. However, if this function is guaranteed to be bounded 
and nowhere zero, then both the multiplication and the division by it are continuous 
operations in either of the function spaces L 2 , C£ and H\ with the understanding that in 
the latter two cases the result belongs to the corresponding spaces if the function itself 
belongs to Cf. In the latter case, by making use of the Neumann series expansion 
the action of the division operator may be traced back to the multiple use of the 
multiplication operator. 

To justify our last claim consider first a continuous operator A acting on a Banach 
space with identity operator /. It is straightforward to check then, by induction, that 

N N 

A - A Y = ~ A Y A = !-(!- A ) N+1 (C- 1 ) 

i=0 i=0 

holds for arbitrary no n- negative integer N. If ||7 — A\\ < 1 we also have that the series 
N ^ Ei=o( J - A Y is absolute convergent (thus its limit is a continuous operator), 
and ||(J — A) Ar+1 || < \\I — A\\ N+1 tends to zero as N — > oo. Thus the relation 
A^ 1 = ^2°^ (I — A) 1 — referred as the Neumann series expansion — follows. Assume now 
that the Banach space in question is C*(M, C), and the operator A is the multiplication 
by a function F G C*(M, C), and denote by ||7 — F\\ the pertinent operator norm of 
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the function 1 — F. Accordingly, an upper bound for the error that we introduce by 
replacing the division operation by multiplication based on the Neumann series may be 



given as 



< 111 



N 

"Ed- 

i=0 

- F\\ N+1 



N 



i - £(i - -fT* 1 



i=0 



(C.2) 



When F is specified by its multipole coefficients, it is not economical to determine the 
norm of 1 — F, as it would require pointwise evaluation of its multipole series. However, 
if F G CI for r > ^ + k, and our Riemann manifold is compact, the Sobolev embedding 
theorem significantly simplifies the determination of an upper bound of the uniform C k 



norm as 



|1 



lc 6 fc 



\1-F\ 



< cm 



(C.3) 



where C is the minimal Sobolev constant, and the Sobolev norm on the right hand side 
may directly be determined by making use of the multipole coefficients of F as 



£ Sr(\) \ (Yi, l) L 2 - (Yi, F) L2 \ 
iex 



(C.4) 



where the function S r introduced in |Appendix A has been applied. Given the value of 
C || 1 — F\\h2 < 1, a predefined error tolerance e can be guaranteed to hold simply by 
calculating Neumann series up to the order N £ = int {ln(e)/ln(C ||1 — F||#2)}, where 
int {x} denotes the integer part of x G K. It is straightforward to see that the number 
of orders to be calculated grows only logarithmically with the increase of the desired 
accuracy. 

The above described method based on the use of the Neumann series may further 
be optimized by rescaling our field F with a complex number 2 in a way to minimize 
(|| 1 — zF\\h2^ . It can be justified, by a direct calculation, that the unique minimum 
may be achieved by choosing 



z(F) 



£, 6X S r (A,)|(^,F> i2 | 
and the corresponding minimal value is 



(C.5) 



M(F) = £5 r (A i )|(K i ,l> 



E, e x Sr(^j) (Yj, F) L2 (Y j: l) 



L 2 



L 2 \ 



E k& xS r (Xk)\(Y k ,F) 

L-'l 



(C.6) 



If we denote the constant eigenfunction of the Laplace operator by Yq, then z(F) and 



M(F) can be re-expressed as z(F) 



1 (Y ,F), 2 



and M(F) 



\\F-{Y ,F) l2 Y P 



£1 j. fly 

optimal rescaling is applied, a sufficient condition for the Neumann series to converge 
is that the inequality C ■ a/ M(F) < 1, where C stands for the Sobolev constant, holds 
and the minimum number of orders necessary to achieve an accuracy below a pre-fixed 
value e is N £ = int{ln(e)/ln(C yJM(F))}. This requirement may also be rephrased 
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as follows. The Neumann series after optimal rescaling is absolute convergent in the 

\\F— (Yo,F) 3 Vq [I 2 

uniform C k norm if rrW^ — — which is nothing but the if ^-measure of the non- 

monopole content in F — is smaller than the threshold As it will be demonstrated by 
the following two examples relevant for the case of the two-spheres and n-dimensional 
toruses, this is a rather weak condition. It is also worth to note that the Neumann 
series expansion may be re-expressed in an iterative form, which requires less function 
evaluations than the canonical series expansion representation does [55J. 



Appendix D. The Sobolev Constants on Two-sphere and on n-Toruses 



As it follows from the discussions in Appendix C to be able to have accurate estimates 



of certain errors, it is also important to know the numerical value of the minimal Sobolev 
constant. In general, the determination of the value of the minimal Sobolev constant 
is a delicate issue. A powerful method yielding this constant is based on the use of 
reproducing kernel property j56] of Sobolev spaces which can be outlined as follows. 
Consider a Hilbert space J4? of some complex valued functions over some set X. Then, 
if the point evaluation / (->■ f(x) is a continuous linear map for every x G X it is called 
a reproducing kernel Hilbert space. The Riesz representation theorem ensures that for 
each x G X there exists a unique K x G such that (K x , f) = f(x) for any / G J4f. As 
K x itself is a function, it may also be evaluated at any point. The reproducing kernel 
function K : X x X i— > C is defined as K(x, y) = K x (y). It may be verified that for any 
x, y G X 

(i) (K(x,-),K(y,-)) = K(y,x), 

(ii) K(y,x) = K(x,y), 

(iii) if comprises a complete orthonormal system, then, K(x, ■) = ^2 ieX <&i(x)<f>i, 
where the infinite summation makes sense in the norm topology. 

If H^(M, C) is a Sobolev space over a compact Riemann manifold M with r > | + k, then 
by the Sobolev theorem the H% norm is stronger than the uniform C k norm, therefore the 
point evaluation is a continuous map. Thus, for an arbitrary choice of an L 2 -orthonormal 
eigensystem {Yi}j g j of the Laplace operator, there exists a unique reproducing kernel 
function K r which — for any choice of x G M reads — as K r (x,-) = Yliex s fy) ^i( x )^ii 
and for any / G H^(M, C) we have (K r (x, •), f)ji 2 (MC) = f( x )- Armed with this identity, 
the relation 

(((V« 9 id) K r ) (x, •), f) HHMC) = (V«/) (x) (D.l) 

can be seen to hold, where the operator (V^ <8> id) acts on K r as the Z-times gradient— 
with 1 = 0,..., k — on the first variable of K r while the second variable of K r remains 
intact. Note that, in virtue of Lebesgue's dominated convergence theorem, the order of 
the action of the gradient V and the scalar product (•, •) may be interchanged. All these 
observations imply that 

|V«/| (x) 
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< 



((V (0 ® id) K r ) I 



x. 



■If), 



'H*{M,C) 

(D.2) 

where in the second step the Cauchy-Schwartz inequality in H^(M,C) has been 
applied. This inequality is known to be sharp, i.e. it may be saturated. When the 



norm || ((V w <g) id) K r ) (x, -)| 



H?(M,C) 



is guaranteed to be independent of x G M — this 



happens, e.g. in case of homogeneous manifolds — we get the sharp inequality 



sup 

|V (0 /| (x) < sup ||((V (Z) ® id) K r ) 



X. 



\H?(M,C) 



H?(M,l 



(D.3) 



Applying this relation to the case of k = the minimal Sobolev constant may be read 
off the particular form of (D.3) as 



c° 



< sup ||i^ r (x, 



\H?(M,C) 



ff*(M,C) 



(D.4) 



where r > | is tacitly assumed to hold. Note, however, that for higher value of k the 
above argument does not necessarily lead to a sharp inequality. 

Let us restrict again considerations to the case of a two- sphere where, according to 
the above discussion, the reproducing kernel can be given as 

^ ' l{t + 1) - 1 



EE 



-Y t m (x)Y r e 



(D.5) 



IS 



Then, as the two- sphere is a homogeneous manifold || ((V w ® id) K r ) (x, •) 
independent of the location of x on S> 2 . By choosing x to be the north pole in standard 
spherical polar coordinates, and also by using the values of the spherical harmonics at 
the north pole we immediately get that 

_1_^(2^ + 1)(^ + 1)-1) 

47T 



\K r (x, 



1=0 



{^ + i)Y+ 1 - 1 



(D.6) 



The square root of the right hand side of (D.6 ) provides the minimal value of the Sobolev 
constant, C r , over the two-sphere with k = and r > 1. The approximate numerical 
values of this Sobolev constants C r for the particular values of r = 2, 3, 4 are listed in 
Table [QU 

Sobolev constant 



i- 



1.284533 



■ 1.106732 
i= • 1.048986 

47T 



Table Dl. Approximate values of the minimal Sobolev constant in the C° C 
Sobolev embedding over the two-sphere, for the r = 2, 3, 4 values. 



Let us finally restrict attention to the case of an n-dimensional torus. Then, the 
reproducing kernel may be given as 

+°° +°° v <t\ V 

Kr {x,)= E ••• E si-v- -£v (D - 7) 
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where ||((V (/) <g> id) K r ) (x, -)\\ H 2 is constant as a function of x G T n , as the n-torus is 
also a homogeneous manifold. By choosing x to be the point where all the polar angle 
coordinates are zero we immediately get 

+00 +OG 



\K r 



(27T) 



E - E 



1 



S r ( k\ ... k„) 



(D.8) 



fei=— 00 k n = 

The square root of this expression gives the minimal Sobolev constants, C r , for the 
Sobolev inequality with r > | and k = 0. For the particular case of T 2 , the approximate 
numerical value of the Sobolev constants C r for k = and for the particular values of 
r = 2, 3, 4 are listed in Table |D2 



Sobolev constant 
J ■ 1.943685 
^ • 1.547391 
f • 1.397749 



Table D2. Approximate values of the minimal Sobolev constant in the C° C 
Sobolev embedding over the two-torus, for the r = 2, 3, 4 values. 



Appendix E. The Estimation of the Tail Sum Error 

As our numerical method is based on the use of multipole expansion of the basic field 
variables, and also since, in practice, we always use only a finite number of multipole 
components, it is of crucial importance to provide precise estimates on the pertinent 
errors. An immediate upper bound on the Sobolev norm of the truncated part of a 
function / may be given as follows. 



Let / G H%(M, C) for some r, then, in virtue of (A. 14), we have that 



^Sr^KXuf)^ =\\f\\ 2 m <00, (E.l) 

iex 

which implies that the sequence % h-> S r (Xi) |(1^,/) L2 | 2 is summable. Consider now a 
sequence of positive numbers i 1— > such that the relation 

\(Y l ,f) L2 \<a l (E.2) 

holds for each i G Z. Then, because the summation preserves monotonicity, the tail sum 
of the sequence i 1— > S r (\i) |a;| 2 bounds the if 2 norm-square of the tail sum error of /. 
Such a bounding sequence may be readily constructed by assuming that / G H ^ (M, C) 
for some r' > r, which implies that also 

Y,sA\)mj) L2 \ 2 = 11/11^ < 00 (e.3) 

holds. Then, the sequence i 1— > dj may be chosen to be an arbitrary monotonically 
decreasing sequence, for which the sum J^ez ^r'(^i) \ a i\ 2 is divergent, as in that case, 
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there always exists a threshold index, above which the relation (E.2) holds. Then 



i i — y o>i may be normalized in such a way that (E.2) holds for any i EX. Such a minimal 
multiplier would be max — \ (Yi, f)\, which, in practice, may always be identified by the 
pertinent maximum on the stored finite orders. Note that this approximation is exact 
whenever the threshold index is reached within the stored orders. 

Restricting again considerations to the case of a two-sphere, a suitable bounding 
sequence (£, m) !->■ a™ may be chosen as 

where K is an unknown normalization factor. Then, the inequality 



2^ 2^ uu + !))-! I^/'^I 

<|K| 2 V ^ £ + 1 )) r+1 - 1 1 fE5) 
(£(£ + l)) r '+ 1 -12£ + l 1 J 

holds for the if 2 norm-square of the tail sum, whenever \K\ is chosen to be large enough 

such that the inequality \(XT^f)i 2 \ — l-^lsS+i ( ^(i^i+i))-]" 1 ) * s satisfied for any 
allowed values of £, m, where r < r'. An immediate lower bound for such \K\ may be 
given as 



£=o77oo m=-"w V (-^(^ + 1)) - 1 



max max LI , ^ — (2€+ 1) |(^, /> i2 | . (E.6) 



By construction, the maximum value in question is necessarily attained at certain finite 
but unspecified I. In numerical applications, an estimate for \K\ may be obtained by 
assuming that the stored components up to the order £ max already encode the relevant 
multipole orders for the field, i.e., the multipole coefficients in the tail already fit into 
the above scheme of the asymptotics. In this case, we may choose K as 



K -JS3L fiyyiTV + 1) I OT. fi„ I ■ (e.7) 

It follows from the above observations that whenever / G C°°(§ 2 ,C) then — since 
for any r f G i? 2 /(§ 2 ,C) also holds — its multipole coefficients (Y™,f) L2 shall decay 
faster than any polynomial in £. Note finally, that this upper bound on the if 2 norm 
of the tail sum error may be used, along with the Sobolev constant, to derive an upper 
bound on the C° norm of the tail sum error. 



Appendix F. The Evaluation of the Integrals determining Radiation 
Anisotropy 

In various physical applications it is important to determine the flux of conserved 
quantities through a hypersurface with topology KxC. 
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In the applied framework of multipole expansions the calculation of these type of 
quantities is rather straightforward, as, in general, they may be written as integrals of 
multilinear expressions of the basic variables over the compact surface C. Thereby they 
can be given in terms of L 2 (C, C) scalar products of some of the multilinear expressions 
of the basic field variables. Recall that multilinear expressions may be evaluated by 



purely spectral methods — as pointed out in Appendix A and Appendix B — , whereas 
the L 2 scalar product may be calculated directly, using the fact that the eigenfunctions 
of the Laplace operator comprise an orthonormal system. 

However, when we are interested in the scale of the anisotropy in the distribution 
of certain quantities, e.g. the energy radiated inwards or outwards, which cannot be 
expressed via the aforementioned integrals over the entire compact surface C. Instead, 
the pertinent integrals have to be evaluated over a subset A C C to compare the radiated 
flux through A and its complement in C. These type of integrals may, however, be given 
as (•, Xa')l 2 (c c) sca l ar products of multilinear expressions of basic field variables over 
C, where \ A denotes the characteristic function of the set A. As the multiplication by 
X A is a continuous map, all that we need for the spectral evaluation of this expression 
is the values of the (V^, X A ^j) l 2 (c c) e m& ^ x elements. 

To demonstrate that this process is much more straightforward in practice than it 
sounds let us restrict our considerations again to the case of the two-sphere. Then we 
have that 



/V m l -v V m 2\ — / V mi V ma 

\ 1 h > X A 1 e 2 I ip- ~ / h 1 i-2 

J A 

min(fl,£2) 

/ i \rr12 \ , /i-mi,m2,-(m 2 -mi) / ™-mi (V \\ 



k=0 " A 

where the expansion of products of spherical harmonics was used, and the symbol G 



stands for the Gaunt coefficients (B.l) 



Let us further restrict considerations to an axially symmetric subset A C C. Then, 
in our parametrization, the integral J A Y^ 2 _^^_ 2k vanishes, unless m 2 = m\. Moreover, 
the integral f A Y^ i _ g ^ +2k may be calculated by making use of the identity 

P' e+1 -Pi 1 = (2£+l)P i . (F.2) 

This relation, which is valid for Legendre polynomials, may easily be verified by using 
the Rodriguez formula. As a consequence, the relation 

/ Y° e = v^o - ./-^— [P m (costf) - P,_x(costf)] (F.3) 

J[0,i?]x[0,27r] V liL + L 

holds, where 5 denotes Kronecker delta. Note that in evaluating the Legendre 
polynomials well-known computational methods can be applied (for a relevant 
implementation see, e.g. [54]). 



References 



[1] Pretorius F 2005 Evolution of binary black hole spacetimes Phys. Rev. Lett. 95 121101 



Multipole Expansion in Time Evolution of Non-linear Dynamical Systems 



47 



[2] Campanelli M, Lousto C 0, Marronetti P and Zlochower Y 2006 Accurate evolutions of orbiting 

black hole binaries without excision Phys. Rev. Lett. 96 111101 
[3] Baker J G, Centrella J, Choi D I, Koppitz M and van Meter J 2006 Gravitational wave extraction 

from an inspiraling configuration of merging black holes Phys. Rev. Lett. 96 111102 
[4] Shibata M and Nakamura T 1995 Evolution of three-dimensional gravitational waves: Harmonic 

slicing case Phys. Rev. D 52 5428-5444 
[5] Baumgarte T W and Shapiro S L 1999 On the numerical integration of Einstein's field equations 

Phys. Rev. D 59 024007 

[6] Loffler F, Faber J, Bentivegna E, Bode T, Diener P, Haas R, Hinder I, Mundim B C, Ott C D, 
Schnetter E, Allen G, Campanelli M, Laguna P 2011 The Einstein toolkit: a community 
computational infrastructure for relativistic astrophysics Preprint gr-qc/1 11 1.3344 
[7] Shibata M and Taniguchi K 2011 Coalescence of black hole-neutron star binaries 

http:/ /www. livingreviews.org/lrr-201 1-6 
[8] Burko L M and Khanna G 2009 Late-time Kerr Tails Revisited Class. Quantum Grav. 26 015014 
[9] Tiglio M, Kidder L and Teukolsky S A 2008 High accuracy simulations of Kerr tails: coordinate 
dependence and higher multipoles Class. Quantum Grav. 25 105022 
[10] Scheel M A, Erickcek A L, Burko L M, Kidder L E, Pfeiffer H P and Teukolsky S A 2004 3D 

simulations of linearized scalar fields in Kerr spacetime Phys. Rev. D 69 104006 
[11] Zenginoglu A and Tiglio M 2009 Spacelike matching to null infinity Phys. Rev. D 80 024044 
[12] Zenginoglu A 2008 Hyperboloidal evolution with the Einstein equations Class. Quantum Grav. 25 
195025 

[13] Zenginoglu A and Kidder L E 2010 Hyperboloidal evolution of test fields in three spatial dimensions 
Phys. Rev. D 81 124010 

[14] Robertshaw O 2011 A spectral MHD code for the non-linear study of magnetic neutron stars PhD 

thesis, University of Southampton 
[15] Racz I and Toth G Z 2011 Numerical investigation of the late-time Kerr tails Class. Quantum 

Grav. 28 195003 

[16] Williams R K 2004 Collimated escaping vortical polar e~e + jets intrinsically produced by rotating 

black holes and Penrose processes The Astrophys. J. 611 952 
[17] Gariel J, MacCallum A H, Marcilhacy G and Santos N O 2010 Kerr geodesies, the Penrose process 

and jet collimation by a black hole Astronomy & Astrophysics 515 A15 
[18] Takami K and Kojima Y 2009 Collimation of a spherical collisionless particles stream in Kerr 

spacetime Class. Class. Quantum Grav. 26 085013 
[19] Csizmadia P, Laszlo A and Racz I 2010 Linear waves on fixed Kerr background and their relevance 

in jet formation Journ. Phys. Conf. Ser. 218 012007 
[20] The GridRipper 3+ld PDE solver homepage 

http : / / www. rmki . kfki.hu / ~ gridripper 
[21] Teukolsky S A 1972 Rotating Black Holes: Separable Wave Equations for Gravitational and 

Electromagnetic Perturbations Phys. Rev. Lett. 29, 1114-1118 
[22] Press W H and Teukolsky S A 1973 Perturbations of a Rotating Black Hole. II. Dynamical Stability 

of the Kerr Metric Astrophys. J. 185, 649-674 
[23] Lehner L, Palenzuela C, Liebling S.L, Thompson C, Hanna C 2011 Intense Electromagnetic 

Outbursts from Collapsing Hypermassive Neutron Stars Preprint gr-qc/1112.2622 
[24] Etienne Z.B, Liu Y.T, Paschalidis V, Shapiro S.L 2012 General relativistic simulations of black 

hole-neutron star mergers: Effects of magnetic fields Phys. Rev. D 85 064029 
[25] Gustafsson B, Kreiss H O and Oliger J 1995 Time dependent problems and difference methods 

Pure and Applied Mathematics (New York: Wiley) 
[26] Wald R M 1984 General relativity (Chicago: University of Chicago Press) 

[27] Chandrasekhar S 1983 The mathematical theory of black holes (Oxford: Oxford University Press) 
[28] Csizmadia P 2006 Testing a new mesh refinement code in the evolution of a spherically symmetric 
Klein-Gordon field Int. J. Mod. Phys. D 15 107 



Multipole Expansion in Time Evolution of Non-linear Dynamical Systems 



48 



Csizmadia P 2007 Fourth order AMR and nonlinear dynamical systems in compactified space 

Class. Quantum Grav. 24 S369 
Berger M J and Oliger J 1984 Adaptive mesh refinement for hyperbolic partial differential equations 

J. Comput. Phys. 53 484 
Friedrich H and Nagy G 1999 The initial boundary value problem for Einstein's vacuum field 

equations Commun. Math. Phys. 201 619-655 
Friedrich H and Rendall A D 2000 The Cauchy problem for the Einstein equations Led. Notes 

Phys. 540 127-224 

Dafermos M and Rodnianski I 2004 A note on boundary value problems for black hole evolution 
Preprint gr-qc/0403034 

Friedrich H 2009 Initial boundary value problems for Einstein's field equations and geometric 

uniqueness Gen. Rel. Grav. 41 1947-1966 
Sarbach O 2007 Absorbing boundary conditions for Einstein's field equations J. Phys. Conf. Ser. 

91 012005 

Allen E W, Buckmiller E, Burko L M and Price R H 2004 Radiation tails and boundary conditions 

for black hole evolutions Phys. Rev. D 70 044038 
Csizmadia P and Racz I 2010 Gravitational collapse and topology change in spherically symmetric 

dynamical systems Class. Quantum Grav. 27 015001 
Carter B 1968 Global structure of the Kerr family of gravitational fields Phys. Rev. 174 1559 
Penrose R 1969 Gravitational collapse: the role of general relativity Riv. Nuovo Cimento 1 special 

number 252-276 

Krivan W, Laguna P, Papadopoulos P and Andersson N 1997 Dynamics of perturbations of rotating 

black holes Phys. Rev. D 56 3395 
Andersson N, Laguna P and Papadopoulos P 1998 Dynamics of Perturbations of Rotating Black 

Holes. II. A note on superradiance Phys. Rev. D 58 087503 
Starobinskii A A 1973 Amplification of waves during reflection from rotating "black hole" 

Zh. Eksp. Teor. Fiz. 64 48-57 
Finster F, Kamran N, Smoller J, Yau S T 2009 A rigorous treatment of energy extraction from 

rotating black hole Commun. Math. Phys. 287 829-847 
Dafermos M and Rodnianski I 2010 The black hole stability problem for linear scalar perturbations 

Preprint gr-qc/1010.5137 
Messiah A 1961 Quantum mechanics I. (Amsterdam: North-Holland Publ. Comp.) 
Genz A and Malik A 1980 An adaptive algorithm for numerical integration over an n-dimensional 

rectangular region J. Comp. Appl. Math. 6 295 
Price R H 1972 Nonspherical perturbations of relativistic gravitational collapse 1. scalar and 

gravitational perturbations Phys. Rev. D 5 2419 
Price R H 1972 Nonspherical perturbations of relativistic gravitational collapse 2. integer-spin, 

zero-rest-mass fields Phys. Rev. D 5 2439 
Adams R A 1975 Sobolev spaces (Academic Press) 

Aubin T 1998 Some nonlinear problems in riemannian geometry (Springer) 

Slater C 1960 Quantum theory of atomic structure (McGraw-Hill) 

Gaunt J A 1929 On the triplets of helium Philos. Trans. Roy. Soc. A 228 151 

Messiah A 1962 Quantum Mechanics (North-Holland) 

The GNU Scientific Library 

|http: / / www.gnu.org/ software/gsl| 
Laszlo A 2006 A robust iterative unfolding method for signal processing J. Phys. A 39 13621 
Aronszajn N 1950 Theory of reproducing kernels Transactions of the American Mathematical 

Society 68 337 



